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CHAPTER  1 


Introduction 


1.1  Dissertation  Objectives 

The  objectives  of  this  research  are  to  further  evaluate  and  improve 
Goddard  Space  Flight  Center’s  (GSFC’s)  Global  Positioning  System  (GPS) 
Enhanced  Orbit  Determination  Experiment  (GEODE)  [1]  software.  GEODE  is 
evaluated  with  various  processing  scheme  changes  with  the  understanding  that  it 
will  be  employed  for  real-time  satellite  Orbit  Determination  (OD),  on  board  a 
satellite.  Improvements  to  GEODE  are  sought  in  terms  of  enhanced  autonomy, 
improved  accuracy/precision  and  reduced  computational  burden.  The  goal  of  this 
research  is  to  autonomously  process  GPS  pseudoranges,  in  real-time,  to  produce 
orbit  estimates  with  better  than  1 -meter  total  position  (3D)  Root  Sum  Square 
(RSS)  error  for  Low  Earth  Orbit  (LEO)  satellites.  The  estimates  produced  by 
GEODE  are  compared  to  Precise  Orbit  Ephemeris  (POE)  generated  by  GSFC  and 
the  Jet  Propulsion  Laboratory  (JPL). 

1.2  Motivation 

The  topic  of  accurate,  autonomous  real-time  satellite  OD  is  receiving 
significant  attention  [2-8]  and  with  GPS  Selective  Availability  (SA)  being  turned 


off  on  2  May  2000,  it  will  receive  even  more  attention  in  the  near  future.  Mission 
planners  want  very  accurate  (tens  of  cm  to  tens  of  meters)  position,  velocity 
and/or  attitude  information  in  real-time,  while  minimizing  the  work  required  to 
achieve  these  results  (better,  faster,  cheaper).  Several  trends  driving  autonomy, 
timing,  and  accuracy  requirements  are  planned  deployment  of  constellations  of 
satellites,  the  desire  for  real-time  geodetic  measurements,  challenging  Earth 
resource  science  objectives  and  a  movement  toward  reducing  dependence  on 
ground-based  tracking  assets  [9].  It  is  even  becoming  desirable  to  perform  OD  in 
real-time,  on  board  an  Earth  orbiting  satellite,  where  accurate  position,  velocity 
and  attitude  information  are  made  available  for  other  satellite  instruments  [4]. 

An  obvious  choice  for  providing  observations  for  an  autonomous  OD 
scheme  is  the  Global  Positioning  System  (GPS).  GPS  provides  unprecedented 
observability  to  LEO  satellites  providing  continuous,  all-weather  observations 
without  user  intervention.  The  limitation  in  using  GPS  for  precise  OD  used  to  be 
the  accuracy  of  the  measurements.  The  User  Equivalent  Range  Error  (UERE)  as 
published  by  the  GPS  Joint  Program  Office  (JPO)  was  33.3  meters  (la)  using  the 
GPS  Standard  Positioning  Service  (SPS)  [10].  SPS  GPS  measurement  errors 
were  dominated  by  errors  due  to  SA  but  SA  has  now  been  turned  off  (as  of  2  May 
2000).  Now  GPS  measurement  errors  are  dominated  by  the  ionosphere, 
troposphere,  GPS  satellite  ephemerides,  multi-path,  and  receiver  noise.  The  effect 
the  measurement  errors  have  on  satellite  OD  has  been  minimized,  if  not 
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completely  removed,  by  post-processing  the  GPS  measurements  using 
sophisticated  filters  and  data  gathered  from  globally  distributed  receiving  stations. 

Therefore,  the  distinction  between  post-processing  (or  processing  in  near 
real-time)  and  real-time  is  very  important  when  discussing  the  accuracy  attained 
using  GPS  measurements  for  satellite  OD.  Post-processing  (also  near  real-time) 
implies  providing  solutions  at  least  several  hours,  usually  more,  after  observations 
are  taken.  Real-time  implies  calculating  solutions  within  minutes,  seconds  or 
even  a  fraction  of  a  second  after  observations  are  made.  Herein,  real-time  OD 
will  be  defined  as  completing  the  calculations  required  to  perform  the  OD 
measurement  update  prior  to  acquiring  the  next  measurement.  The  main 
differences  between  these  two  types  of  systems  are  that  post-processing  systems 
usually  corrected  for  SA,  use  very  sophisticated  ionospheric  models,  use  high 
fidelity  gravity  models  and  use  very  accurate  GPS  satellite  ephemerides  rather 
than  those  broadcast.  Each  of  these  “differences”  require  a  prohibitive 
computational  burden  for  a  real-time  system  or  unacceptable  interaction  for  true 
autonomy.  Obviously,  since  SA  has  been  turned  off,  post-processing  systems  no 
longer  correct  for  it.  Autonomy  defined  here,  implies  minimal  or  no  interaction 
with  external  resources,  except  the  acquisition  of  measurements  and  the  collection 
of  information  in  the  GPS  navigation  message.  The  next  two  sections  of  this 
chapter  outline  several  state-of-the-art  post-processing  and  real-time  software 
suites  to  distinguish  the  differences  between  these  two  types  of  systems. 


1.3  Post-processing  Systems 


There  are  several  software  packages  designed  to  calculate  precise  satellite 
ephemeris  using  GPS  observations  by  post-processing  (on  the  ground)  or  in  near 
real-time.  Several  examples  are  the  Jet  Propulsion  Laboratory’s  (JPL)  GPS- 
Inferred  Positioning  System  (GIPSY)/OASIS  II  (GOA  II)  [11],  the  Naval 
Research  Laboratory’s  (NRL)  Orbit/Covariance  Estimation  and  Analysis 
Software  (OCEAN)  [12],  Van  Martin  System  Inc.’s  (VMSI)  MicroCosm®, 
GEODYN,  Utopia,  etc.  GOA  II,  OCEAN  and  MicroCosm®  are  described  here. 


1.3.1  JPL’s  GOA  II 

GOA  II  is  a  collection  of  Fortran  programs  and  C-shell  or  Perl 
encapsulating  scripts  created  by  JPL.  Most  of  the  following  features  of  GOA  II 
were  extracted  from  [13]: 

•  batch-sequential  and  filter-smoother 

•  network  processing  of  globally  distributed  GPS  receiver  data  to  eliminate 
SA  using  white  noise  clock  estimation  with  one  reference  clock  [14] 

•  simultaneous  estimation  of  GPS  satellite  orbits  (generation  of  precise  GPS 
Ephemeris  (GPSE)) 

•  orbit  integration  using  spherical  harmonic  gravity  field  expansion  (JGM-2 
70x70  or  JGM-3  70x70),  the  effects  of  the  sun,  moon,  and  planets,  plus 
non-gravitational  forces  to  account  for  atmospheric  drag  and  solar 
radiation  pressure 

•  modeling  of  the  known  dynamics  of  the  Earth,  including  solid  Earth  tides, 
precession,  nutation,  polar  motion,  pole  tides,  and  ocean  loading 

•  ionospheric  group  range  delay  and  phase  advance  estimation  based  on 
Bent  or  IRI95  ionosphere  models 

•  wet  and  dry  tropospheric  delay  modeling  including  ray  bending  and  Earth 
curvature,  and  the  stochastic  estimation  of  mapped  zenith  delays  and  clock 
biases  as  random  walk  or  Gauss-Markov  process  noise 


•  estimation  of  unmodeled  or  mismodeled  forces  using  the  Reduced 

Dynamic  Technique  (RDT)  [15] 

GPS  satellite  orbits  computed  using  GOA  II  have  a  radial  RMS  in  the  7-10 
cm  range.  Comparison  of  The  Ocean  Topography  Experiment 
(TOPEX)/Poseidon  (T/P)  orbits  derived  from  GPS  data  collected  on  board  the  T/P 
satellite  compared  with  orbits  produced  by  other  sources  (e.g.,  Doris  and  Satellite 
Laser  Ranging  (SLR)  analyzed  by  other  institutions)  have  shown  2-3  cm  radial 
RMS  and  15  cm  3D  RSS  [13].  JPL  has  determined  operational  GPS  orbits  for 
T/P  for  many  years  with  typical  latency  of  1 1-17  hours  after  the  last  GPS  data 
point  is  received  on  board  the  satellite.  In  comparison,  GSFC  definitive  orbits 
based  on  Doris  and  SLR  data  are  delivered  with  a  latency  of  about  40  days.  The 
definitive  GSFC  orbits  and  JPL  orbits  generated  with  GOA  II  have  an  RMS 
difference  of  less  than  4  cm  radially  [16].  JPL  has  automated  GOA  II  to  produce 
orbits  with  sub-10  cm  radial  RMS  accuracy  with  about  10-hour  latency  but  recent 
improvements  in  JPL’s  processing  system  might  soon  reduce  this  to  around 
1-hour  latency.  JPL  anticipates  OD  accuracy  at  the  1  cm  level  (la  radially)  will 
be  demonstrated  in  the  not  too  distant  future  [16].  In  this  dissertation  GEODE 
T/P  position  estimates  will  be  compared  against  JPL  GOA  II  POE. 

1.3.2  NRL’s  OCEAN 

OCEAN  is  a  precise  satellite  OD  software  system  created  by  the 
Astrodynamics  and  Space  Applications  Office  of  the  Naval  Center  for  Space 


Technology  at  the  NRL.  The  OCEAN  suite  contains  a  batch  weighted  least 
squares  estimation  technique  and  an  extended  Kalman  filter  and  backwards 
smoother  to  process  various  range,  Doppler  and  angle  observations.  It  can 
process  GPS  data  or  SLR  data.  The  previous  information  and  the  following 
features  were  extracted  from  [12]. 

•  Choice  of  Earth  gravity  models  including  JGM-2  or  EGM96 

•  Jaccia  1971  atmospheric  drag  model 

•  Solar  radiation  pressure  using  vehicle  macro-modeling 

•  Sun,  Moon  and  planetary  third  body  accelerations 

•  Models  solid  Earth  and  Ocean  tide  accelerations 

•  Processes  externally  available  precise  or  broadcast  GPS  ephemerides 

•  Capable  of  processing  multiple  satellites  and  ground  stations 
simultaneously 

•  Differential  correction  using  ground  station  receiver  data 

•  Estimation  of  unmodeled  or  mismodeled  (empirical)  accelerations 


OCEAN  was  used  to  process  T/P  data  from  18  November  1993  when  all  GPS 
satellites  had  Anti-Spoofing  (AS)  off  but  only  17  of  24  satellites  had  SA  off.  AS 
is  the  encryption  of  the  GPS  P-code  to  keep  adversaries  from  “spoofing”  a  Precise 
Positioning  System  (PPS)  capable  receiver  [17].  Turning  AS  off  essentially 
allows  a  dual  frequency  receiver,  without  decryption  capabilities,  to  receive 
carrier  phase  and  pseudorange  measurements  on  both  the  LI  and  L2  channels. 
Only  the  data  from  the  17  S  A  off  satellites  were  processed.  The  results  compared 
to  JPL’s  T/P  POE  with  precise  GPSE  are  shown  below  in  Table  1.1. 


Table  1.1  -  OCEAN  OD  Results 


RMS 

Precise  GPSE 

(simulated  post-processing) 

Radial 

0.10  m 

Along-track 

0.31m 

Cross-track 

0.16  m 

Total 

0.37  m 

Again,  the  GPS  data  in  this  study  are  SA  and  AS  free. 


1.3.3  VMSI’s  MicroCosm® 

MicroCosm®  is  a  satellite  orbit  and  geodetic  parameter  estimation 
software  suite  developed  by  VMSI.  MicroCosm®  improves  upon  and  fully 
implements  the  NASA  GEODYN II,  version  8609,  precision  orbit  and  geodetic 
parameter  determination  software  system.  The  following  details  concerning 
MicroCosm®  were  extracted  from  [18]: 


•  Bayesian  least  squares  batch  processor 

•  automatic  carrier  phase  cycle  slip  detection  and  removal 

•  network  processing  of  globally  distributed  GPS  receiver  data  to  eliminate 
SA  using  single,  double  or  triple  differencing 

•  simultaneous  estimation  of  geodetic  parameters 

•  orbit  integration  using  spherical  harmonic  gravity  field  expansion  (JGM-2 
70x70  or  JGM-3  70x70),  the  effects  of  the  sun,  moon,  and  planets,  plus 
non-gravitational  forces  to  account  for  atmospheric  drag  and  solar 
radiation  pressure 

•  modeling  of  the  known  dynamics  of  the  Earth,  including  solid  Earth  tides, 
precession,  nutation,  polar  motion,  pole  tides,  and  ocean  loading 

•  estimation  of  zenith  tropospheric  parameters  for  each  IGS  ground  station 

•  estimation  of  radial,  in-track  and  cross-track  empirical  accelerations 


MicroCosm®  was  used  to  perform  orbit  determination  for  the 
GPS/Meteorology  (GPS/MET)  experiment  using  GPS  carrier  phase  observations. 
Internal  and  external  orbit  overlap  comparisons  show  MicroCosm®  is  capable  of 
estimating  GPS/MET  orbits  at  the  30  cm  3D  RSS  level  [19].  In  this  dissertation, 
GEODE  GPS/MET  position  estimates  will  be  compared  against  MicroCosm® 
generated  GPS/MET  orbits. 

1.4  Real-time  Systems 

There  are  also  several  real-time  software  suites  available.  Examples  are 
JPL’s  Real-time  GIPSY,  the  Microcosm  Autonomous  Navigation  System 
(MANS)  [20],  the  Brazilian  National  Institute  of  Space  Research’s  (INPE’s) 
ORBesT  [21]  and  GSFC’s  GEODE  [22].  Also,  the  University  of  Nottingham’s 
Institute  of  Engineering  Surveying  &  Space  Geodesy  (EESSG)  developed  another 
unnamed  system  [2].  There  is  currently  no  published  information  on  space 
qualified  (actually  flown  in  space)  precise,  real-time  OD  software.  Real-time 
GIPSY  and  the  University  of  Nottingham’s  system  are  discussed  below;  GEODE 
is  discussed  in  Chapter  2. 

1.4.1  JPL’s  Real-time  GIPSY  (RTG) 

RTG  is  an  ANSI  C  version  of  GOA-II  created  by  JPL  to  accommodate 
high  data  rates  (1  Hz)  and  improve  portability  to  systems  other  than  UNIX.  JPL’s 
goal  is  to  incorporate  all  the  precise  models  from  GOA-II,  make  it  suitable  for 


imbedded  systems  such  as  GPS  receivers  and  make  it  capable  of  real-time 
processing  [16].  Compiler  options  in  RTG  allow  it  to  be  scaled  to  meet  various 
processor  load  requirements  [16].  To  provide  the  best  accuracy,  RTG  will  be 
used  in  conjunction  with  a  global  Wide  Area  Augmentation  System  (WAAS)  or  a 
Wide  Area  Differential  GPS  (WADGPS)  system.  Without  WAAS  or  WADGPS 
RTG  has  shown  3D  RSS  values  in  the  4-6  meter  range  when  used  to  process  T/P 
data  with  broadcast  GPS  ephemeris  and  S A  on  [16]. 

1.4.2  The  University  of  Nottingham’s  Study 

In  a  study  for  the  UK  Defense  and  Evaluation  Research  Agency  (DERA), 
The  Institute  of  Engineering  Surveying  &  Space  Geodesy  (IESSG)  at  the 
University  of  Nottingham  developed  an  extended  Kalman  filter,  using  Reduced 
Dynamic  Tracking  (RDT),  to  generate  real-time  satellite  position  estimates  with 
radial  RMS  error  of  1.08  m  1  a  and  a  3D  RSS  error  of  3.95  m  la  [2].  IESSG 
used  real  and  simulated  Standard  Positioning  System  (SPS)  data  from  T/P.  They 
reported  the  filter  converged  after  approximately  five  hours  [23].  They  used  a 
JGM-2  45x45  gravity  field,  a  simple  drag  model  (due  to  T/P’s  relatively  high 
orbit),  and  broadcast  GPS  ephemerides.  The  application  required  approximately 
500  kb  of  computer  memory  and  the  code  could  produce  solutions  within  one 
minute  of  recording  an  observation.  A  trade  study  between  microprocessors  was 
also  performed  finding  a  military  standard  1750A  microprocessor  (8086 
equivalent)  to  be  more  than  capable  of  producing  the  solutions  each  minute  [2]. 
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1.5  Post-processing  /  Real-time  Comparison 

A  comparison  of  the  characteristics  common  to  post-processing  and  real¬ 
time  systems  is  shown  below  in  Table  1.2.  Details  of  research  outlined  in  this 
dissertation  to  close  the  gap  between  post-processing  and  real-time  systems  are 
presented  next. 


Table  1.2  -  Post-processing/Real-time  Comparison 


Area  of  Concern 

Post-processing 

Real-time 

SA  Treatment 

Differential  Correction 

Broadcast  DGPS 

GPSE  Used 

Precise 

Broadcast 

Gravity  Model 

JGM-2  70x70,  JGM-3 
70x70  or  EGM-96 
300x300 

Truncated  Models  (45x45 
largest,  usually  30x30) 

Drag  Model 

Complex 

Simple 

Radiation  Pressure  Model 

Complex 

Simple  or  Not  at  All 

Ionospheric  Model 

Bent,  IRI95,  PRISM 

None 

Observations  Used 

Pseudorange,  Integrated 
Doppler 

Pseudorange 

Filter  Type 

Usually  Batch  or 
Sequential 

Extended  Kalman  Filter 

Empirical  Force  Estimation 

Usually  RDT 

None 

1.5.1  Error  Sources  in  Real-time 

SA  was  the  U.S.  Air  Force’s  intentional  degradation  of  GPS  measurement 
accuracy.  SA  was  formally  implemented  on  25  March  1990  and  turned  off  on  2 
May  2000.  SA  is  accomplished  through  the  dithering  of  GPS  satellite  clocks  [10] 
The  usual  method  for  dealing  with  SA  is  to  collect  data  at  geographically 
separated  locations  with  known  position  and  use  single  or  double  differencing  to 
remove  SA  errors.  Fortunately,  SA  is  no  longer  a  factor  in  real-time  satellite  OD 
as  there  is  currently  no  real-time  method  to  reduce  its  contribution  to  degrading 


OD  accuracy. 


Another  error  source  that  is  difficult  to  deal  with  in  real-time  is  error  in  the 


knowledge  of  the  GPS  satellite’s  position.  Currently  the  method  used  to  calculate 
GPS  satellite  position,  in  real-time,  is  to  evaluate  a  series  of  equations  using 
information  broadcast  by  each  GPS  satellite.  In  Zumberge  [24]  the  RMS  errors 
between  the  broadcast  and  precise  GPS  ephemeris  for  the  period  1  July  1993  -  22 
Oct  1993  are  reported  as  1.2  m  radial,  3.2  m  cross-track  and  4.5  m  along-track. 
The  total  3D  RSS  error  is  therefore  5.7  m.  Depending  on  the  type  of  GPS  data 
utilized  (including  or  not  including  SA  and/or  AS),  the  difference  in  accuracy 
between  the  broadcast  and  precise  ephemerides  has  been  reported  to  be  a 
significant  source  of  error  in  determining  the  orbit  of  a  satellite  using  GPS 
observations  [12, 25-27].  To  mitigate  the  errors  in  broadcast  GPS  ephemerides, 
precise  ephemerides  are  often  used.  There  are  two  methods  of  obtaining  the 
precise  ephemerides:  solving  for  the  GPS  orbits  simultaneously  with  the  user 
satellite  state  and  fixing  the  GPS  ephemerides  to  an  independent  determination, 
such  as  the  International  GPS  Service  for  Geodynamics  (IGS)  [28].  Both  of  these 
methods  require  significant  overhead  and  are  therefore  not  implemented  in  real¬ 
time  systems.  However,  predicted  GPS  ephemerides  are  available  and  in  certain 
circumstances  could  be  uplinked  daily  to  a  satellite.  Chapter  2  presents  a 
comparison  of  real-time  satellite  OD  error  results  using  broadcast  and  precise 
GPS  ephemerides. 

The  model  used  to  predict  the  acceleration  due  to  the  Earth’s  gravity  can 
be  a  significant  error  source  in  real-time  systems.  Therefore,  a  tradeoff  study  is 


performed  in  Chapter  3  to  determine  the  effect  of  using  various  truncated  gravity 
models.  In  addition  a  gravity  acceleration  approximation  method  is  introduced  to 
recover  the  accuracy  of  large  models  while  reducing  computational  burden. 

Measurement  errors  due  to  the  ionosphere  can  also  affect  satellite  OD 
accuracy.  Chapter  4  presents  details  concerning  the  magnitude  of  errors  due  to 
ionospheric  path  delay  seen  by  the  GPS/MET  satellite  and  an  innovative  method 
of  removing  these  errors  when  a  single  frequency  GPS  receiver  is  used.  The 
technique  presented  is  known  as  Differenced  Range  Verses  Integrated  Doppler 
(DRVID)  [29], 

Due  to  limited  space  environment  knowledge,  computational  burden 
restrictions  and  the  desire  for  autonomy,  all  the  forces  acting  on  a  satellite  cannot 
be  modeled.  Therefore,  Dynamic  Model  Compensation  (DMC)  is  presented  in 
Chapter  5  to  not  only  estimate  unmodeled  and  mismodeled  forces  as  part  of  the 
filter  state  but  also  to  simultaneously  provide  a  method  of  formulating  the  filter 
process  noise  variance-covariance. 

One  challenge  with  DMC  is  tuning  of  the  time  correlation  coefficient  and 
the  process  noise  standard  deviation.  Chapter  6  presents  a  Genetic  Algorithm 
(GA)  in  the  form  of  Genetic  Model  Compensation  (GMC)  [9].  GMC  adaptively 


tunes  the  constants  used  in  DMC. 
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1.6  Dissertation  Summary 

There  are  significant  differences  between  the  satellite  OD  accuracy 
attainable  in  a  post-processing  scenario  as  compared  to  in  real-time.  The  purpose 
of  the  research  presented  in  this  dissertation  is  to  offer  suggested  improvements  to 
GSFC’s  GEODE  software  suite  to  enhance  autonomy,  improve  accuracy  and 
reduce  computational  burden. 

Chapter  1  provides  a  definition  of  the  term  real-time  and  a  review  of  state- 
of-the-art  post-processing  and  real-time  software  suites  and  their  differences. 
Chapter  2  discusses  details  of  GEODE  and  the  spacecraft  and  data  used  to 
evaluate  GEODE’s  performance.  Chapter  3  deals  with  the  accuracy  and 
computational  burden  of  the  gravitational  model  used  in  GEODE  and  suggests  a 
promising  method  to  significantly  reduce  the  computational  burden.  This 
gravitational  approximation  sacrifices  computer  memory  but  not  accuracy. 
Chapter  4  describes  the  errors  due  to  the  ionosphere  in  GPS  measurements  and 
suggests  a  method  to  estimate  and  remove  ionospheric  error  using  a  technique 
called  Differenced  Range  Versus  Integrated  Doppler  (DRVID)  [29-32]. 

Chapter  5  presents  a  derivation  of  an  XYZ  and  a  Radial,  In-track,  Cross-track 
(RIC)  version  of  Dynamic  Model  Compensation  (DMC)  [9],  results  of  two 
simulations  implementing  RIC  DMC  and  the  improvements  to  GEODE  through 
the  use  of  DMC.  Chapter  6  details  an  extension  of  DMC  by  the  use  of  a  Genetic 
Algorithm  (GA)  to  adaptively  tune  the  time  correlation  coefficient  (t  )  and 


standard  deviation  of  the  DMC  forcing  noise  ( a )  [9].  Finally,  Chapter  7  provides 
conclusions  and  suggestions  for  future  work. 


CHAPTER 2 


GEODE  and  GPS/MET 

This  chapter  provides  an  overview  of  the  OrbView-1  satellite  (formerly 
known  as  MicroLab-1)  carrying  the  GPS/MET  experiment  and  information 
concerning  the  GPS  measurements  collected  by  GPS/MET.  Details  are  provided 
on  GEODE  and  results  presented  concerning  the  optimization  of  inputs  to 
GEODE.  Finally,  a  Lagrange  scheme  is  introduced  to  interpolate  precise  GPS 
ephemerides  and  a  method  presented  for  removing  the  effects  of  SA  from 
GPS/MET  data  collected  2-10  Feb  1997. 

2.1  GPS/MET 

OrbView-1  was  launched  on  a  standard  Pegasus  rocket  from  Vandenberg 
AFB,  CA,  on  3  April  1995  into  a  747  x  732  km  orbit  with  an  inclination  of  70.0°. 
The  primary  payload  on  board  OrbView-1  is  the  Optical  Transient  Detector 
(OTD)  used  in  imaging  lightning  strikes.  OTD  has  a  10-km  resolution,  a 
1,300-km  swath  width,  and  is  managed  by  the  NASA  Marshall  Space  Flight 
Center  [33].  GPS/MET  is  a  secondary  payload  on  board  the  OrbView-1  satellite 
built  by  Orbital  Sciences  Corporation  (OSC).  The  GPS/MET  experiment  relies 
on  an  active  limb  sounding  technique  using  radio  occultation  observations  taken 


by  the  on  board  GPS  receiver.  Figure  2.1  shows  a  diagram  of  the  occultation  of  a 
GPS  signal  by  the  Earth’s  atmosphere. 

Non-Occulting 

Occulting  GPS 
Satellite 

Atmosphere 


Figure  2.1  -  Occultation  Description 

The  GPS/MET  receiver  is  a  JPL  modified  TurboRogue  with  a  microstrip 
patch  antenna  [19].  The  occultation  data  is  used  for  recovery  of  accurate 
refractivity,  pressure,  temperature  and  moisture  profiles  of  the  atmosphere  [34]. 
Optimal  occultation  tracking  is  performed  when  the  GPS/MET  antenna  is 
pointing  in  the  anti-velocity  direction,  therefore,  a  fixed  yaw  steering 
configuration  is  implemented.  OrbView-1  is  gravity-gradient  stabilized  and 
attitude  is  maintained  using  three  torque  rods,  six  Sun  sensors,  two  Earth  sensors 
and  a  magnetometer  mounted  on  the  gravity  boom.  Attitude  accuracy  of  5-10°  is 


maintained.  The  GPS/MET  antenna  phase  center  is  offset  from  the  OrbView-1 
satellite’s  center  of  gravity  (CG)  by 


-0.507  m  in  the  x  direction 
0.039  m  in  the  y  direction 
-0.178  m  in  the  z  direction  [19] 


See  Figure  2.2  for  a  diagram  and  Figure  2.3  for  a  photo  of  OrbView-1. 


Figure  2.2  -  OrbView-1  and  GPS/MET 


/ 


Figure  2.3  -  OrbView-1  Satellite 

Photo  Courtesy  of:  http://www.orbimage.com/satellite/orbviewl/orbviewl.html 


2.2  GEODE 


GSFC  developed  GEODE  as  a  real-time  software  analysis  package  [35]. 
GEODE  is  highly  modular,  programmed  in  ANSI  C  and  has  been  targeted  to 
UNIX  and  PC  systems  as  well  as  the  RAD6000  RISC  microprocessor.  It  requires 
a  modest  400  kilobytes  of  computer  RAM.  GEODE  was  originally  designed  as 
experimental  software  to  fly  on  the  SSTI  Lewis  satellite  contracted  by  NASA  to 
TRW  [35].  GEODE  is  implemented  with  an  Extended  Kalman  Filter  (EKF), 
which  feeds  a  real-time  state  propagator.  GEODE  is  designed  to  be  hosted  on 
either  a  spacecraft  flight  computer,  or  in  a  GPS  receiver's  processing  unit.  Pre¬ 
launch  orbit  determination  studies  using  GEODE  indicate  that  1  a  orbit  accuracy 
of  10  m  in  position,  and  0.01  m/s  in  velocity  may  be  attained  in  the  presence  of 
SA.  The  pre-launch  studies  were  accomplished  with  Extreme  Ultraviolet 
Explorer  (EUVE)  and  T/P  raw  pseudorange  data.  The  SSTI  Lewis  satellite  was  to 
be  placed  in  a  523  km.  Sun  Synchronous  orbit  [36].  Unfortunately,  SSTI  Lewis 
was  lost  shortly  after  its  launch  in  August  of  1997  and  therefore  GEODE  is  not 
flight  qualified.  Below  is  a  summary  of  relevant  information  concerning  GEODE 
[37], 

•  JGM-2  30x30  gravity  model 

•  Solar  and  Lunar  point  mass  3rd  body  force  model 

•  Harris-Priester  atmospheric  drag  model 

•  Geometrical  editing  of  measurements  with  high  ionospheric  errors 

•  Broadcast  GPS  ephemerides  used 

•  Extended  Kalman  Filter  (EKF)  implemented 

•  UDUt  factorized  state  error  covariance 


•  Uplink  of  polar  motion  coefficients,  accurate  a  priori  state  information  and 
a  priori  state  error  and  process  noise  covariance  terms 

•  Processes  pseudorange  measurements  only  -  not  carrier  phase 

Figure  2.4  is  generated  using  position  data  estimated  by  GEODE 
compared  to  post  processed  POE  (accurate  to  approximately  30  cm  radially) 
generated  with  MicroCosm®.  The  observations  used  by  GEODE  are  GPS 
pseudoranges  collected  by  the  GPS/MET  experiment  on  4  Feb  1997.  AS  is  off 
but  SA  is  on.  There  are  24  hours  of  data  with  one  pseudorange  processed  every 
10  seconds.  It  took  less  than  one  minute  to  process  the  entire  24  hours  of  data  on 
a  450  MHz  Pentium  II  with  128  MB  RAM.  The  filter  converged  after  processing 
approximately  two  hours  of  data  and  yielded  a  converged  3D  RSS  error  of 
11.61  m.  All  RMS  and  RSS  statistics  presented  in  this  dissertation  are  calculated 
using  converged  estimates  (after  2  hours)  only.  Figure  2.4  shows  plots  of  the 
Radial,  In-track  and  Cross-track  (RIC)  errors  compared  to  the  MicroCosm®  mean 
of  J2000  POE,  along  with  the  square  root  of  the  filter’s  estimated  state  error 
covariance  (positive  and  negative),  also  called  the  estimated  standard  deviation. 
Figure  2.4  also  shows  the  RIC  RMS  values  and  the  mean  of  the  errors  in  each 
direction.  A  brief  description  of  why  the  RMS  and  mean  of  errors  are  shown  in 
this  dissertation  as  well  as  a  brief  discussion  of  systematic  errors  follows. 

Without  any  systematic  errors  the  RMS  values  presented  would  indicate 
the  level  of  precision  and  accuracy  attained  by  GEODE.  In  the  presence  of 
systematic  errors,  the  RMS  really  only  provides  a  measure  of  precision.  In  Figure 
2.4  systematic  errors  can  be  seen  in  the  periodic  nature  of  the  plots  and  possibly  in 


the  mean  of  the  comparison  error.  For  instance,  the  cross-track  plot  in  Figure  2.4 
shows  a  periodicity  with  frequency  of  once  per  satellite  orbit.  This  “once  per  rev” 
periodicity  is  probably  caused  by  GEODE/POE  coordinate  transformation 
differences  but  could  also  be  due  to  errors  in  the  filter’s  dynamic  and 
measurement  models.  Periodicities  in  errors  are  not  revealed  in  the  mean  and 
RMS  of  error  statistics.  Looking  at  plots  and/or  plotting  the  power  spectrum  of 
the  errors  can  reveal  periodic  errors.  The  in-track  plot  in  Figure  2.4  shows  a  mean 
in  the  difference  between  the  POE  and  the  filter’s  estimate  of  -1.12  m  in  the  in¬ 
track  direction.  It  is  difficult  to  ascertain  the  cause  of  the  error  mean  but  it  could 
be  from  systematic  or  random  errors.  Systematic  errors  can  be  introduced  into  the 
filtering  process  through  the  GPS  pseudorange  measurements,  the  mechanics  of 
the  filter,  the  observation  model  or  the  dynamic  model.  Random  errors  are  most 
certainly  introduced  in  the  GPS  pseudorange  measurements  but  could  also  be 
introduced  in  other  ways. 

Since  the  source  of  the  mean  of  comparison  error  shown  in  the  in-track 
plot  in  Figure  2.4  cannot  be  determined,  no  conclusions  can  be  drawn  from  it 
alone.  However,  if  a  change  is  made  to  GEODE  that  decreases  the  error  mean, 
then  it  might  be  concluded  the  change  reduced  some  systematic  error.  In  this 
dissertation  improvements  to  GEODE  will  be  gauged  not  only  by  their  affect  on 
the  RMS  and  RSS  of  the  errors  between  GEODE’s  satellite  position  estimates  and 
POE  but  also  by  the  affect  on  the  mean  of  the  comparison  errors. 
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Figure  2.4  -  GEODE  GPS/MET  Position  Error  in  Mean  of  J2000  Rotated  to  RIC 
Statistics  Are  Calculated  Based  On  Converges  Estimates  Only 


2.2.1  UDUt  Factorized  Covariance 

The  conventional  Extended  Kalman  Filter  (EKF)  formulation  is  shown 
below  [38]. 


Pk  —  ®(lk>lk-l)Pk-l®  (tkAk-l)  +  Qk 


yk=Yk 

^  observed  •*' 


computed 


Kk  =pI[h^(h(.pi;h^  +r)  1 


+ -^kYk 

Pk=(l-KkHk)Pk 


(2.1) 

(2.2) 

(2.3) 

(2.4) 


(2.5) 


where 


Pk  =  time  updated  covariance  matrix  at  time  k 
Pk_i  =  measurement  updated  covariance  at  time  k-1 
Qk  =  process  noise  covariance  matrix  at  time  k 
yk  =  measurement  residual  at  time  k 
Kk  =  Kalman  gain  at  time  k 
Hk  =  observation-state  matrix  at  time  k 
Xk  =  A  priori  estimate  of  the  state  at  time  k 
Xk  =  best  estimate  of  the  state  at  time  k 
R  =  Observation  Error  Covariance  Matrix 

Implementation  of  the  conventional  EKF  can  lead  to  filter  divergence  either 
through  inaccuracies  in  the  mathematical  models  used  in  the  filter  (dynamic  or 
measurement  models)  or  through  the  state  error  covariance  matrix  becoming  non¬ 
positive  definite.  Since  errors  introduced  during  the  computational  procedure  are 
the  cause  of  the  covariance  matrix  becoming  non-positive  definite,  a 
reformulation  of  the  algorithm  can  minimize  these  errors  [38].  Therefore,  several 
filter  algorithm  modifications  have  been  developed  to  improve  the  numerical 
stability  of  the  Kalman  filter. 

The  UDUt  factorized  covariance  implementation  of  the  Kalman  filter  was  first 
introduced  by  Thornton  and  Bierman  [39, 40]  where 

P  =  UDUt  (2 

U  is  unit  upper  triangular  and  D  is  diagonal.  Factorizing  the  state  error 
covariance  matrix  in  this  way  avoids  square  roots,  guarantees  non-negativity  of 
the  computed  covariance  and  keeps  the  filter  numerically  stable  and  accurate  [37], 


Numerical  stability  means  that  an  algorithm  computes  the  same  result  even  when 
the  initial  conditions  are  slightly  perturbed.  The  UDUT  factorization  provides 
numerical  stability  [41].  The  UDUT  factorization  algorithm  performs  the  time 
and  measurement  updates  on  the  U  and  D  matrices  rather  than  directly  updating 
the  P  matrix.  The  algorithm  is  implemented  in  GEODE  so  the  measurement 
update  is  performed  on  each  measurement  independently.  Maybeck  shows  both 
the  numerical  advantages  and  computational  burden  changes  due  to  employing 
various  covariance  factorizations.  Maybeck’ s  conclusion  is  that  with  a  10 
member  state  and  two  measurements  to  update,  UDUT  factorization  provides  the 
best  balance  of  numerical  stability  and  reduced  computational  burden  [41].  The 
algorithm  as  implemented  in  GEODE  can  be  found  in  the  GEODE  Mathematical 
Specification  [37]. 

2.2.2  Process  Noise 

Kalman  filter  algorithms  use  the  process  noise  covariance  matrix,  Qk,  to 
correct  the  state  error  covariance  for  inadequacies  in  the  force  model.  See 
equation  (2.1).  Without  process  noise,  elements  of  the  covariance  matrix 
asymptotically  approach  zero,  as  does  the  Kalman  gain.  When  this  occurs,  as 
seen  in  equation  (2.4),  the  filter  begins  ignoring  measurements  [38].  Ignoring 
measurements  forces  the  filter  to  rely  solely  on  the  dynamic  model  and  since  the 
system’s  dynamics  are  linearized,  filter  divergence  can  result. 


GEODE’s  EKF  uses  a  “physically  connected”  algorithm  for  calculating 
the  gravitational  acceleration’s  contribution  to  the  position  and  velocity 
components  of  the  state  error  covariance  [37].  The  “physically  connected” 
algorithm  uses  gravity  model  variances  in  the  position  and  velocity  process  noise 
formulations.  Gravitational  acceleration  forces  are  chosen  since  they  provide  the 
largest  force  model  errors.  The  formulation  of  the  position  and  velocity  process 
noise  components  (upper  left  6x6  of  the  Q  matrix)  rely  heavily  on  radial,  in-track, 
and  cross-track  correlation  times  provided  to  GEODE  in  the  uplink  command  file. 
These  correlation  times  change  depending  on  the  semi-major  axis  of  the  user 
satellite’s  orbit  and  the  degree  and  order  of  the  gravity  model  used  in  GEODE. 
The  correlation  times  are  calculated  using  an  algorithm  presented  in  Wright  [42] 
and  Fortran  code  is  available  with  GEODE  for  these  calculations.  Example 
correlation  times  for  various  truncations  of  the  JGM-2  and  EGM-96  gravity 
models  are  shown  in  Chapter  3. 

Random-walk  algorithms  are  implemented  for  the  remaining  elements  of 
the  state,  i.e.,  the  receiver  clock  bias,  receiver  clock  bias  drift  rate,  atmospheric 
drag  coefficient  (Cd)  and  the  solar  radiation  pressure  coefficient  (Cr).  The 
process  noise  covariance  algorithms  for  these  elements  of  the  state  are  also  found 
in  the  GEODE  Mathematical  Specification  [37]. 
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2.2.3  Dual  Frequency  GPS  Observations 

GEODE  was  designed  for  use  with  a  single  frequency  GPS  receiver. 
Therefore,  several  changes  are  made  to  GEODE’s  measurement  handling 
procedures  to  take  advantage  of  dual  frequency  GPS  observations.  The  reason  for 
the  change  is  that  the  GPS/MET  data  was  collected  using  a  dual  frequency 
receiver,  with  AS  off. 


2.2.3.1  Ionospheric  Correction 

If  dual  frequency  GPS  measurements  are  used  in  GEODE,  the  PI  and  P2 
measurements  can  be  combined  to  remove  the  first  order  ionospheric  effects  on 
the  pseudorange  and  phase  measurements  [10, 43].  The  ionospheric  time  delay  at 
the  Li  frequency  is 


At,  — 


1 


f  #2  \ 


A-i 2J 

Likewise  the  ionospheric  time  delay  at  L2  is 


(Ppi  ) 


AtL2=~ 

2  c 


1  f  *2 


TZ7T  (Pr.-PrJ 
1  L2  J 


(2.7) 


(2.8) 


where 


c  =  the  speed  of  light  in  a  vacuum  =  299,792,458  m/s 
fi  =  1575.42  MHz 
f2  =  1227.6  MHz 

pPi  =  pseudorange  measurement  taken  on  Li 
Pp2  =  pseudorange  measurement  taken  on  L2 
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Therefore,  the  ionosphere  free  pseudorange  measurement  is: 

Pp3  =  Pp,  +  =  Pp2  c^l2 


or 


leading  to 


Pp3  —  Pp, 


f r  -  f : 


(pp.-ppj 


2  J 


Pp3  “ 


f  g 


/  f2  A 

12 

f2  —  f2 

VT1  r2  ) 


Pp2 


(2.9) 


(2.10) 


(2.11) 


Substituting  in  the  LI  and  L2  frequencies  and  simplifying  we  get: 


Pp3  - 


1698pPi  -1031pp2 
667  ~ 


an  approximation  to  this  equation  is: 

Pp3  =  2-5pPj  —  1.5pP2 


(2.12) 


(2.13) 


Likewise,  the  following  can  form  the  ionospheric  free  phase  measurement: 


As  ApPi  =  cAtLi , 


leading  to 


and  thus 


A<t>L,  =-fLlAtLl  and  Ac()L2  =-fL2AtLj 

<t>L3I  =  4>l,  and  <j)L32  =  (|)L2  -  Mu_ 

<t>L31  =‘t>L,  +fL,AtLi  and  <t)L3,  =  ^L,  +fL,AtLl 


(2.14) 


(2.15) 


(2.16) 


using  equation  (2.7) 
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Figure  2.5  shows  a  plot  of  the  dual  frequency  ionospheric  corrections  for 
the  4  Feb  1997  GPS/MET  data.  The  solar  flux  (F10.7)  for  the  4  Feb  1997  is 
approximately  70.  The  benign  solar  flux  environment  contributes  significantly  to 
the  low  (3.93  m)  mean  of  the  ionospheric  delay  on  this  day.  The  elevation  of  the 
GPS  satellites  with  respect  to  the  OrbView-1  satellite  is  also  shown  on  the  plot  in 
Figure  2.5  for  the  largest  ionospheric  corrections.  Distinction  between 
measurements  taken  while  the  OrbView-1  satellite  is  in  the  Sun  (day)  and  in 
darkness  (night)  is  also  made  on  the  plot.  The  difference  in  the  mean  between  the 
day  and  night  corrections  is  significant  (2.79  m),  as  expected. 


28 


mean  3.93  m 


E 

•£" 

o 

o 

£ 

o 

O 


CD 

3 

cr 

<D 


CO 

3 

Q 


Figure  2.5  -  GPS/MET  Dual  Frequency  Ionospheric  Delay  Corrections 

2.2.3.2  Interfrequency  Bias 

Another  important  consideration  when  using  dual  frequency  GPS  data  is 
to  correctly  handle  the  GPS  satellite  interfrequency  bias.  Differences  in  the 
hardware  signal  paths  of  the  LI  and  L2  signals  on  board  each  GPS  satellite  create 
a  space  vehicle  (SV)  dependent  group  delay  differential  between  LI  and  L2  [44]. 
The  L1-L2  interfrequency  bias  is  broadcast  in  the  navigation  message  and  is 
designated  with  the  variable  Tg(i.  The  Tgd  correction  is  implemented  for  single 
frequency  users  because  the  broadcast  clock  correction  coefficients  are  based  on 
the  effective  code  phase  using  dual  frequency  ionospheric  corrections.  Therefore, 
the  single  frequency  user  must  adjust  for  this  differential  delay.  GEODE  does 
account  for  Tgd  so  it  must  be  removed  when  processing  dual  frequency  data. 


Table  2.1  shows  a  comparison  of  the  errors  between  GEODE  and  MicroCosm® 
POE  when  GEODE  uses  single  and  dual  frequency  measurements.  Although  the 
3D  RSS  and  mean  of  the  radial  error  errors  improved  with  use  of  the  dual 
frequency  ionospheric  corrections,  the  mean  of  the  in-track  error  changed  by  over 
3  meters.  One  possible  explanation  is  an  unknown  systematic  error  is  introduced 
by  the  application  of  the  dual  frequency  correction.  There  could  be  errors 
associated  with  the  Tgd  correction  itself  as  JPL  found  errors  in  the  initial  broadcast 
Tgd  values  and  began  estimating  new  values  for  the  Air  Force  to  broadcast  starting 
in  1999  [44].  Another  possibly  error  source  could  be  that  the  receiver 
interfrequency  bias  is  not  being  accounted  for  adequately  in  the  estimation  of  the 
receiver  clock  bias.  No  additional  insight  into  the  cause  of  the  in-track  error  mean 
increase  is  found  by  comparing  the  measurement  residuals  in  Table  2.2  or  by 
looking  at  plots  of  the  errors  in  Figure  2.6  and  Figure  2.7.  Qualitatively,  it  can  be 
said  that  the  precision  of  the  GEODE  estimate  using  dual  frequency 
measurements  is  improved  but  nothing  definitively  can  be  said  about  its  accuracy. 
The  dual  frequency  GPS/MET  ionospheric  corrections  will  be  used  in  Chapter  4 
to  compare  GEODE  results  achieved  with  application  of  the  Differenced  Range 
Versus  Integrated  Doppler  (DRVID). 


Table  2.1  -  GEODE  GPS/MET  Single  Versus  Dual  Frequency  Error  Results 


GPS/MET -4  Feb  1997 

Error  Mean  (m) 

RMS  Error  (m) 

R 

I 

c  , 

R 

I 

c 

3D 

Single  Frequency 

0.55 

-0.99 

-0.10 

2.74 

10.24 

2.15 

10.81 

Dual  Frequency 

0.11 

2.04 

-0.11 

2.48 

9.68 

2.12 

10.21 

Cross-track  Error  (m)  In-track  Error  (m)  Radial  Error  (m) 


Table  2.2  -  GEODE  GPS/MET  Single  Versus  Dual  Frequency  Measurement 

Residuals 


GPS/MET  4  Feb  1997 

Measurement 
Residual  RMS  (m) 

Measurement 
Residual  Mean  (m) 

Single  Frequency 

44.96 

-0.11 

Dual  Frequency 

44.60 

-0.07 

3D  RSS  Error  =  10.81  m 


2  4  6  8  10  12  14  16  18  20  22  24 


2  4  6  8  10  12  14  16  18  20  22  24 


4  8  12  16  20  24 

Time  (hours) 


Figure  2.6  -  GPS/MET  Single  Frequency  Results 
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3D  RSS  Error  =  10.21  m 


2  4  6  8  10  12  14  IS  18  20  22  24 


Figure  2.7  -  GPS/MET  Dual  Frequency  Results 


2.2.4  Height  Of  Ray  Path  (HORP)  Editing 

GEODE  uses  a  technique  called  HORP  to  edit  out  measurements  that 
travel  through  the  atmosphere  [37].  HORP  editing  requires  two  user  inputs 
contained  in  the  uplink  command  file:  atmosphere  height,  h,  and  maximum 

central  angel,  otmax.  h  and  ow  are  shown  in  Figure  2.8. 
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Measurements  are  edited  out  if 


d|  <  RE  +  h  and  a  >  .  Here 


a  =  cos 


-i 


|  *user  1 1  %] 


IPS!  J 


and 


elevation  =  —  -  cos  1 
2 


Ax-r  A 

A  \iser 


m 


VI  I I  user |  J 


If  the  input  for  h  in  the  uplink  command  file  is  set  to  0.0  then 


h=kser|~15km 

If  h  is  set  to  0.0  when  processing  GPS/MET  data  all  measurements  taken  below 
approximately  -3.72°  elevation  will  be  considered  for  editing.  Again, 


(2.21) 


(2.22) 


(2.23) 


measurements  will  not  be  edited  out  unless  d  <  RE  +  h  and  a  >  .  As 


shown  in  Table  2.3,  the  value  of  ainax  can  significantly  affect  GEODE’s  results. 


Table  2.3  -  GEODE  GPS/MET  Varying  amax  Results 


GPS/MET -4  Feb  1997 

Mean  Error  (m) 

RMS  Error  (m) 

C*max*  (dCg) 

elevation 

R 

I 

C  . 

R 

I  ■;  . 

c 

3D 

70 

4.6 

0.12 

4.11 

-0.07 

4.18 

14.77 

2.64 

15.58 

75 

-0.5 

0.12 

4.11 

0.07 

4.18 

14.77 

2.64 

15.58 

80 

-5.4 

0.09 

3.49 

-0.08 

4.46 

15.41 

2.64 

16.26 

85 

-10.2 

0.04 

2.71 

-0.14 

4.15 

13.35 

2.62 

14.23 

90 

-14.9 

0.24 

0.12 

-0.13 

3.81 

11.71 

2.30 

12.52 

95 

-19.6 

0.43 

-1.90 

-0.13 

3.15 

11.17 

2.25 

11.82 

96 

-20.5 

0.42 

-1.54 

-0.11 

3.11 

11.04 

1.99 

11.64 

97 

-21.4 

0.48 

-1.06 

-0.11 

2.86 

10.63 

1.92 

11.18 

98 

-22.3 

0.54 

-1.12 

-0.12 

2.94 

10.96 

2.42 

11.61 

99 

-23.2 

0.63 

-1.34 

-0.11 

3.09 

11.05 

2.22 

11.69 

100 

-24.1 

0.55 

-0.99 

-0.10 

2.74 

10.24 

2.15 

10.81 

101 

-25.0 

0.67 

-1.94 

-0.09 

2.96 

10.53 

2.01 

11.12 

102 

-25.9 

0.66 

-1.61 

-0.09 

3.08 

10.53 

2.01 

11.15 

103 

-26.8 

0.96 

-3.18 

-0.09 

3.01 

11.06 

1.86 

11.61 

104 

-27.7 

0.99 

-3.71 

-0.09 

2.99 

11.27 

2.07 

11.85 

105 

-28,5 

0.99 

-3.71 

-0.09 

2.99 

11.27 

2.07 

11.85 

*  is  nominally  set  to  70° 


**  h  =  0.0  implies  15km  below  the  satellite  altitude  of  ~747  km 


2.2.5  Polar  Motion  and  AUT1 

In  determining  the  transformation  matrices  between  the  Earth  Centered 
Inertial  (ECI)  mean  of  J2000  coordinate  system  and  an  Earth  Centered  Earth 
Fixed  (ECEF)  coordinate  system,  two  parameters  are  used  which  require 
prediction,  and  therefore,  are  challenging  to  implement  in  a  real-time  filter  flown 
on  board  a  satellite.  The  two  parameters  are  the  pole  location  (polar  motion)  and 
the  difference  between  the  UT1  and  UTC  time  systems. 
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2.2.5.1  Polar  Motion 

Polar  motion  takes  into  account  the  fact  that  the  celestial  ephemeris  pole, 
which  is  normal  to  the  true  equator,  is  in  motion  with  respect  to  the  terrestrial 
reference  frame.  In  other  words,  the  rotational  pole  moves  within  the  Earth.  The 
maximum  amplitude  of  polar  motion  is  0.3  arc  seconds  or  approximately  9  m  on 
the  surface  of  the  Earth  [45].  Polar  motion  consists  largely  of  two  motions,  an 
annual  elliptical  component  with  a  period  of  365  days  and  a  Chandler  circular 
component  with  a  period  of  about  435  days  [46].  The  motion  is  difficult  to 
predict  and  is  determined  by  observations. 

The  coordinates  of  the  Earth’s  instantaneous  pole  location  are  measured 
by  the  International  Polar  Motion  Service  (IPMS)  in  terms  of  xp  and  yp 
components  in  the  polar  plane.  xp  is  measured  along  the  Greenwich  meridian  and 

yp  is  measured  along  the  90°  W  meridian.  Past  values  of  xp  and  yp  are  published 

in  the  International  Earth  Rotation  Services  Final  Bulletin  found  at: 
ftp://maia.usno.navv.mil/ser7/finals.all 

A  plot  of  xp  and  yp  from  May  1976  to  May  2000  are  shown  in  Figure  2.9. 
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Figure  2.9  -  Polar  Motion  Plot 

The  United  States  Naval  Observatory  (USNO)  also  supplies  polar  motion 
predictions  published  in  the  International  Earth  Rotation  Service  (IERS)  Bulletin- 
A.  The  estimated  accuracy  of  the  predictions  is  0.005  arc-second  for  the  40-day 
predictions  [37].  The  predicted  values  of  xp  and  yp  are  obtained  by  evaluating  the 
following  trigonometric  functions  for  the  day  of  interest: 

xp  =  a7  +  a2  cos  A  +  a3  sin  A  +  a4  cos  B  +  a5  sin  C  (arc-seconds)  (2.24) 

yp  =  a6  +  a7  cos  A  +  a8  sin  A  +  a9  cos  B  +  a10  sin  C  (arc-seconds)  (2.25) 

where 

A  =  3^25  (MJD  “ Tp  )  (radians) 


(2.26) 


C  =  ^L(MJD-Tp)  (radians) 
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(2.27) 


MJD  =  Modified  Julian  Date  =  JD  of  interest  -  2500000.5  days 
Tp  =  epoch  of  the  prediction 


The  10  coefficients  are  published  in  the  IERS  Bulletin-A  found  at: 


http://hpiers.obspm.fr/iers/bul/bula/bulletinA 


The  effect  of  not  including  polar  motion  in  GEODE  depends  on  the  size  of 
the  pole  wander.  The  results  in  Table  2.4  show  the  effects  on  GEODE’s  position 


estimates  with  and  without  polar  motion  when  GPS/MET  data  is  processed.  The 


comparison  is  again  made  against  one  day  of  MicroCosm®  POE. 


Table  2.4  -  GEODE  GPS/MET  Results  With  and  Without  Polar  Motion 


GPS/MET -4  Feb  1997 

Mean  Error  (m) 

RMS  Error  (m) 

R 

I 

C 

R 

V  I:-:,,,:-; 

C 

3D 

With  Polar  Motion 

0.55 

-0.99 

-0.10 

2.74 

10.24 

2.15 

10.81 

Without  Polar  Motion 

0.65 

-0.99 

-0.12 

2.98 

13.05 

3.16 

13.76 

2.2.S.2  AUT1 

The  calculation  of  the  Greenwich  Hour  Angle  (GHA)  requires  Greenwich 
Sidereal  Time  (GST),  which  uses  the  UT1  time  system.  UT1  is  a  time  scale 
measured  by  the  rotation  of  the  Earth.  UTC  is  the  time  scale  used  worldwide  for 
technical  and  scientific  activities  and  is  a  compromise  between  highly  stable 
atomic  time  and  irregular  Earth  rotation.  The  current  practice  is  to  keep  the 
difference  between  UT1  and  UTC  less  than  0.9  seconds  by  adjusting  UTC  by 


integer  leap  seconds.  The  USNO  distributes  the  UT1-UTC  prediction  known  as 
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AUT1  in  the  IERS  Bulletin-A.  AUT1  is  calculated  by  the  following  equation  as 
presented  in  the  IERS  Bulletin-A: 

AUT1  =  u1+u2(MJD-Tuti)  +  u3(MJD-Tut1)2  (seconds)  (2.28) 


where 


MJD  =  Modified  Julian  Date  =  JD  of  interest  -  2500000.5  days 
Tun  =  MJD  of  the  epoch  of  prediction 


The  estimated  accuracy  of  the  AUT1  prediction  is  0.0048  seconds  for  a  40  day 
prediction.  The  effect  of  not  including  AUT1  in  GEODE  is  far  more  significant 
than  not  including  polar  motion.  See  Table  2.4  and  Table  2.5. 


Table  2.5  -  GEODE  GPS/MET  Results  With  and  Without  AUT1 


GPS/MET  -  4  Feb  1997 

Mean  Error  (m) 

RMS  Error  (m) 

R 

I 

c 

R 

■  'I 

C 

3D 

With  AUT1 

0.55 

-0.99 

-0.10 

2.74 

10.24 

2.15 

10.81 

Without  AUT1 

1.68 

25.80 

-0.9 

4.37 

28.96 

59.73 

66.53 

2.2.6  All  in  View  Versus  Cyclic 

GEODE  has  the  capability  to  process  measurements  from  all  GPS 
satellites  in  view  or  to  process  the  measurements  from  only  one  satellite  at  a  given 
observation  epoch.  The  keywords  associated  with  these  different  processing 
schemes  in  the  uplink  command  file  are  ALL  for  all  in  view  and  CYCLIC  for 
processing  only  one  observation  at  each  observation  epoch.  The  reason  the 


CYCLIC  option  is  implemented  is  SA.  Since  the  effects  of  SA  are  time 
correlated  it  is  thought  that  by  spacing  the  processing  of  GPS  observations  at  the 
same  period  as  the  correlation  in  SA,  the  SA  errors  would  appear  more  random 
[4].  Therefore,  a  time  variable  is  also  included  in  the  uplink  command  file  known 
as  minimum  sampling  frequency.  This  minimum  sampling  frequency  is 
nominally  set  at  the  correlation  time  of  SA  (3-5  minutes)  [37].  In  reality,  tuning 
the  minimum  sampling  frequency  yields  significant  changes  in  GEODE’s 
estimates.  Table  2.6  shows  a  comparison  of  GEODE  results  when  processing  all 
in  view  verses  cyclic.  Table  2.7  shows  results  of  tuning  the  minimum  sampling 
frequency.  An  additional  advantage  to  only  processing  one  GPS  measurement  at 
each  observation  epoch  is  reduced  computational  burden  since  measurements 
from  a  minimum  of  5  and  a  maximum  of  8  GPS  satellites  are  available  at  each 
observation  epoch.  The  CYCLIC  scheme  is  used  throughout  this  dissertation  as 
no  accuracy  /  precision  improvement  is  shown  when  processing  all  GPS  satellites 
in  view,  even  with  SA  free  data. 


Table  2.6  -  GEODE  GPS/MET  ALL  vs.  CYCLIC  Comparison 


GPS/MET -4  Feb  1997 

Mean  Error  (m) 

RMS  Error  (m) 

R 

I 

C 

R 

■  i 

C 

3D 

CYCLIC  -  *110 

0.55 

-0.99 

-0.10 

2.74 

10.24 

2.15 

10.81 

ALL 

1.68 

-4.48 

-0.08 

4.38 

13.94 

2.05 

14.76 

*  1 10  is  the  minimum  sampling  frequency  in  seconds 


Table  2.7  -  GEODE  GPS/MET  Min  Sampling  Frequency  Tuning  Results 


GPS/MET -4  Feb  1997 

Mean  Error  (m) 

RMS  Error  (m) 

min  sampling  frequency 

R 

i 

c 

R 

i 

C 

3D 

10  (sec) 

0.67 

-2.26 

-0.06 

3.36 

11.53 

1.77 

12.14 

20 

0.61 

-1.87 

-0.09 

2.99 

10.77 

1.98 

11.35 

30 

0.60 

-1.70 

-0.10 

3.03 

10.83 

2.17 

11.45 

40 

0.60 

-1.70 

-0.10 

3.03 

10.83 

2.17 

11.45 

50 

0.62 

-2.27 

-0.09 

3.01 

10.75 

1.93 

11.33 

60 

0.61 

-1.82 

-0.09 

2.93 

10.67 

2.00 

11.24 

70 

0.66 

-1.84 

-0.09 

3.07 

10.98 

2.07 

11.59 

80 

0.61 

-1.52 

-0.10 

2.70 

10.44 

2.03 

10.97 

90 

0.65 

-1.67 

-0.10 

2.85 

10.43 

2.15 

11.03 

100 

0.62 

-2.04 

-0.10 

3.09 

11.06 

2.25 

11.70 

110 

0.55 

-0.99 

-0.10 

2.74 

10.24 

2.15 

10.81 

120 

0.68 

-1.88 

-0.08 

3.22 

10.83 

1.77 

11.43 

130 

0.63 

-1.88 

-0.10 

3.16 

10.46 

2.31 

11.17 

180 

0.49 

-0.98 

-0.09 

3.16 

10.10 

2.76 

10.94 

240 

0.50 

i  -0.45 

-0.06 

3.98 

12.18 

2.80 

13.12 

300 

0.27 

-0.02 

-0.12 

4.48 

15.09 

3.06 

16.04 

2.2.7  Antenna  Phase  Center  Offset  From  Satellite  Center  of  Gravity  (CG) 

Since  the  phase  center  of  a  GPS  antenna  may  not  be  coincident  with  the 

CG  of  the  satellite  it  is  important  to  account  for  this  offset  when  computing  the 

predicted  measurement.  As  reported  earlier,  the  phase  center  of  the  GPS/MET 

antenna  is  offset  from  the  OrbView-1  CG  by 

-0.507  m  in  the  x  direction 
0.039  m  in  the  y  direction 
-0.178  m  in  the  z  direction  [19] 

The  OrbView-1  satellite’s  attitude  is  maintained  with  5 -10° accuracy.  Since  a 
10°  attitude  error  translates  to  a  maximum  of  a  4  cm  error  in  x,  y,  or  z,  the  offsets 
detailed  above  can  be  assumed  to  be  fixed.  Here,  x  is  in  the  in-track  direction,  y 
is  in  the  cross-track  direction  and  z  is  in  the  radial  direction.  GEODE  currently 


has  the  capability  to  map  the  antenna  phase  center  to  the  satellite’s  CG  through  a 


transformation  from  RIC  coordinates  to  inertial  coordinates.  Obviously,  a  more 


rigorous  approach  would  be  to  use  accurate  attitude  information  along  with  the 


above  body  fixed  offsets  to  more  precisely  account  for  the  offset  in  inertial  space. 


Unfortunately,  GEODE  does  not  currently  have  this  capability.  Results  are 


shown  in  Table  2.8  using  no  offset  and  the  offset  reported  above.  The  26  cm  3D 


RSS  and  over  1  m  in-track  bias  improvements  came  with  no  additional 


computational  burden  as  GEODE  accounts  for  the  offset  even  when  it  is  zero. 


Table  2.8  -  GEODE  GPS/MET  Antenna  Phase  Center  Mapping  Results 


GPS/MET  -  4  Feb  1997 

Mean  Error  (m) 

RMS  Error  (m) 

Antenna  Offset  (m) 

:  R 

I 

... 

R 

I 

c 

3D 

0.0 

0.0 

0.0 

0.63 

-2.06 

-0.10 

2.70 

10.54 

2.03 

11.07 

0.178 

-0.507 

-0.039 

0.55 

-0.99 

-0.10 

2.74 

10.24 

2.15 

10.81 

2.2.8  Lagrange  Interpolation  of  Precise  GPS  Ephemerides 

The  motivation  for  investigating  the  accuracy  of  precise  GPS  ephemerides 
compared  to  broadcast  ephemerides  is  to  determine  the  accuracy  gained  in 
satellite  OD  and  the  feasibility  of  using  precise  GPS  ephemerides  in  a  real-time 
OD  scheme  on  board  a  satellite.  Several  organizations,  JPL  being  one  of  them, 
produce  GPS  ephemeris  predictions  that  are  more  accurate  than  those  broadcast. 
Satellite  OD  in  real-time  can  use  either  broadcast  or  predicted  precise 
ephemerides.  Obviously,  post-processed  precise  ephemerides  cannot  be  used  in  a 
real-time  scenario.  This  section  includes  a  comparison  of  broadcast  and  precise 
ephemerides,  a  comparison  of  GEODE  interpolated  ephemerides  and  ephemerides 


interpolated  using  an  International  GPS  Service  (IGS)  executable  (reported 
accuracy  of  this  interpolator  is  at  the  mm  level),  a  comparison  of  Lagrange 
interpolated  ephemerides  and  IGS  ephemerides  and  finally  the  Lagrange 
interpolation  algorithm. 

2.2.8.1  Broadcast  Versus  Precise  GPS  Ephemerides 

The  first  comparison  is  between  broadcast  and  precise  ephemerides  for 
GPS  PRN01  on  2  Feb  1997.  Zumberge  [24]  reports  broadcast  ephemerides  are 
accurate  to  5  to  10  meters.  Figure  2.10  shows  a  slightly  smaller  3D  RSS  position 
error  of  3.71  m,  however,  only  4  hours  of  data  are  plotted.  The  time  bias  error  in 
the  lower  of  the  two  plots  in  Figure  2.10,  is  the  difference  between  the  broadcast 
GPS  clock  correction  and  the  correction  supplied  with  the  precise  ephemerides. 
The  large  step  decrease  in  the  bias  error  at  the  two-hour  point  is  due  to  switching 
to  a  new  set  of  broadcast  ephemerides. 


Position  RSS  Error  =  3.7072  m 


Figure  2.10  -  Broadcast  vs.  Precise  Ephemeris  PRN01  2  Febl997 

2.2.8.2  GEODE  Interpolation  Versus  Precise  GPS  Ephemerides 

The  next  comparison  is  between  GEODE’s  interpolation  of  precise 
ephemerides  and  interpolation  done  by  the  IGS  software.  “Truth”  data  is 
generated  using  the  IGS  software  and  output  at  1 -second  intervals.  GEODE’s 
interpolator  is  then  used  to  interpolate  between  5-min,  15-min  and  30-min  interval 
data,  simulating  the  uplink  of  a  set  of  precise  ephemerides.  The  30-min  interval 
data  would  be  the  best  choice  since  it  requires  the  least  uplink  overhead.  Figure 
2.11  shows  the  results  of  GEODE’s  interpolator  using  precise  ephemerides  at  a 
5-min  interval.  It  is  not  known  what  causes  the  excursion  at  the  start  of  the 
position  error  plot.  Figure  2.12  shows  the  results  of  GEODE’s  interpolator  with 
15-minute  interval  precise  ephemerides.  GEODE’s  interpolator  is  not 


documented  in  the  Mathematical  Specification  so  it  is  difficult  to  ascertain  the 
cause  of  the  excursions  and  larger  error.  The  15-minute  interval  GPS 
ephemerides  could  not  be  used  due  to  the  large  errors. 


Position  RSS  Error  =  3.8995  m 
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Figure  2.11  -  GEODE  Interpolator  with  5-Minute  Interval  Ephemeris 


44 


Position  RSS  Error  =  629.0342  m 


Bias  RMS  Error  =  0.0006  m 


Figure  2.12  -  GEODE  Interpolator  with  15-Minute  Interval  Ephemeris 


2.2.8.3  Lagrange  Interpolation  Versus  Precise  GPS  Ephemerides 

Next,  a  Lagrange  interpolation  scheme  is  used  and  the  results  compared  to 
the  IGS  interpolation.  The  best  results  are  attained  using  the  5-minute  ephemeris 
interval  but  the  result  achieved  with  the  12th  order  interpolation  and  30-minute 
interval  is  the  best  balance  of  accuracy,  upload  minimization  and  reduced 
computational  burden.  It  is  interesting  to  note  that  the  Lagrange  interpolation 
performed  worse  using  17th  order  than  12th  order.  Figure  2.13  shows  the  12th 
order  Lagrange  interpolation  using  30-minute  interval  data  and  Figure  2.14  shows 


Position  Error  (m) 


the  17th  order  interpolation  using  30-minute  interval  data.  A  summary  of  the 
results  for  each  of  the  above  methods  is  shown  in  Table  2.9. 
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0.006 

0.004 

0.002 
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Figure  2.13  -  Lagrange  with  30-Minute  Interval  and  12th  Order  Interpolation 


0.015 


Position  RSS  Error  =  0.0055  m 
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Bias  RMS  Error  =  0.00343  m 


Figure  2.14  -  Lagrange  with  30-Minute  Interval  and  17th  Order  Interpolation 


Table  2.9  -  Statistics  of  Results 


Method 

Ephemeris 
Interval  (min) 

Interpolation 

Order 

RSS  of  Range 
Error  (meters) 

RMS  of  Bias 
Error  (meters) 

Broadcast 

n/a 

n/a 

3.7072 

0.05129 

GEODE 

5 

n/a 

3.8995 

0.0001556 

GEODE 

15 

n/a 

629.03 

0.0006 

Lagrange 

30 

12 

0.0040 

0.001692 

Lagrange 

30 

17 

0.0055 

0.00343 

2.2.8.4  Precise  Ephemerides  Interpolation  Conclusions  and  Lagrange 
Algorithm 

The  conclusion  drawn  from  these  results  is  that  the  GEODE  interpolator 
needs  to  be  replaced.  It  produces  inconsistent  results  and  it  requires  ephemerides 
closely  spaced  in  time.  The  Lagrange  interpolation  scheme  performs  more 
consistently,  is  much  more  accurate  with  ephemerides  at  lower  frequency  and  its 


computational  burden  is  smaller  than  GEODE’s.  Below  find  the  algorithm  used 
for  Lagrange  interpolation  borrowed  from  Hofmann-Wellenhof  [17]. 
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Assume  functional  values  f(tj)  are  given  at  epochs  tj,  j  =  0,  1, n.  Where 
n  is  the  order  of  the  interpolation.  Here  f(tj)  are  the  GPS  satellite  x,  y,  z  position 
values  and  the  GPS  clock  bias  values.  Then, 


i  ft)  (j  t»)(t  V) 

The  interpolated  value  at  epoch  t  follows  from  the  summation 

j=0 

The  implementation  of  the  algorithm  shown  above  in  GEODE  is: 

Given:  GPS  ephemeris  (Xj,  yj,  Zj)  and  bias;  at  times  tj  i  =  1  to  n 
Find:  GPS  ephemeris  (x,  y,  z)  and  bias  at  time  T 

x  =  0  y  =  0  z  =  0  bias  =  0 

i  =  1  to  n 

li  =  1 

h  =  time  of  ephemeris  at  time  i 
j  =  1  to  n 

ifi 

tj  =  time  of  ephemeris  at  time  j 
li  =  li  *  (T  -  tj)  /  (ti  -  tj) 

end  if 

endj 

X  =  X  +  li  *  X; 

y  =  y  +  lj  *  yj 

Z  =  Z  +  lj  *  Zj 

bias  =  bias  +  1;  *  bias* 


(2.29) 


(2.30) 


end  i 
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2.2.8.5  GEODE  Results  Using  Precise  GPS  Ephemerides 


There  is  only  a  marginal  improvement  in  GEODE’s  position  estimates, 
compared  to  MicroCosm®  POE,  when  precise  GPS  ephemerides  are  used  instead 
of  those  broadcast.  Table  2.10  shows  a  comparison  of  the  results  using 
GPS/MET  data.  The  minimal  improvement  is  most  likely  due  to  SA.  SA  can 


cause  tens  of  meters  of  error  in  the  pseudorange  measurement  and  is  such  a 
dominating  error  source  it  swamps  the  errors  introduced  by  using  broadcast 
ephemerides. 


Table  2.10  -  GEODE  GPS/MET  Results  Using  Precise  GPS  Ephemeris 


GPS/MET -4  Feb  1997 

RMS  Error  (ra) 

R 

i 

c 

R 

i 

C 

3D 

Broadcast  Ephemeris 

0.55 

-0.99 

-0.10 

2.74 

10.24 

2.15 

10.81 

Precise  Ephemeris 

0.58 

-1.05 

-0.09 

2.99 

10.11 

1.97 

10.73 

Table  2.1 1  shows  results  of  processing  GEODE  on  SA  free  T/P  data  from 
5  May  2000.  Here  the  GEODE  solution  is  compared  to  POE  generated  by  JPL. 
While  the  improvement  on  the  GPS/MET  data  is  only  0.74%,  the  improvement  to 
the  T/P  data  is  10.0%.  Therefore,  it  is  concluded  that  in  the  presence  of  SA  using 
precise  GPS  ephemerides  is  unnecessary  while  in  the  absence  of  SA,  using 
precise  GPS  ephemerides  may  be  a  viable  method  of  improving  OD  accuracy. 


Table  2.11  -  GEODE  T/P  Results  Using  Precise  GPS  Ephemeris 


T/P -5  May  2000 

RMS  Error  (m) 

R 

i 

C 

R 

i 

c 

3D 

Broadcast  Ephemeris 

0.21 

-0.04 

0.01 

0.35 

1.06 

0.44 

1.20 

Precise  Ephemeris 

0.22 

-0.16 

0.00 

0.34 

0.97 

0.31 

1.08 
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2.2.9  JPL  High-Rate  GPS  Clock  Estimates 

In  order  to  assess  GEODE’s  accuracy  with  SA  off,  data  with  SA  free 
measurements  are  needed.  The  T/P  data  mentioned  in  the  previous  section  can  be 
used  but  at  the  T/P  altitude  (1336  km),  ionospheric  effects  are  smaller  and  the 
satellite’s  dynamic  environment  is  more  benign.  GPS/MET,  on  the  other  hand,  is 
lower  (750  km  altitude)  and  thus  experiences  larger  measurement  errors  due  to 
ionosphere  and  a  more  challenging  dynamic  environment.  Unfortunately,  there 
are  no  periods  during  the  GPS/MET  mission  when  SA  is  off.  Therefore,  high  rate 
GPS  clock  estimates  produced  by  JPL  using  GOA  II  are  used  instead  of  the 
broadcast  or  precise  GPS  clock  estimates  [47].  The  high  rate  clock  estimates  are 
found  at: 

ftp://sideshow.ipl.nasa.gov/pub/gipsv  products/hrclocks/ 

2.2.9.1  JPL  High-Rate  GPS  Clock  Estimate  Challenges 

There  are  several  challenges  in  using  the  JPL  high  rate  GPS  clock 
estimates.  There  are  sporadic  time  periods,  for  different  GPS  satellites,  when  the 
high  rate  clock  estimates  are  not  available.  Figure  2.15  shows  two  several  hour- 
long  periods  where  the  high  rate  estimates  are  not  available.  The  3  Feb  1997 
GPS/MET  data  has  only  7  out  of  25  GPS  satellites  without  gaps.  The  typical  gap 
length  is  2-3  hours. 
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Another  challenge  of  using  the  high  rate  clock  estimates  is  there  are 
discontinuities  from  one  day  to  the  next  due  to  changes  in  the  reference  ground 
clocks  used  to  generate  the  estimates  [47].  Also,  satellite  clock  estimates 
available  on  one  day  may  not  be  available  the  next  day.  Fortunately,  each 
24-hour  arc  uses  the  same  reference  ground  clock  and  usually  satellites  are  not 
dropped  during  a  24-hour  arc. 

A  final  item  of  interest  is  a  bias  between  the  broadcast  clocks  and  the  high 
rate  clocks.  As  seen  in  Figure  2.16  and  Figure  2.17,  the  high  rate  GPS  clock 
estimates  from  3-5  Feb  1997  are  biased  by  72.3  m  for  PRN  2  and  71.8  m  for 
PRN  9.  The  bias  is  not  constant  over  the  coarse  of  a  day  and  slowly  increases 


over  time.  The  data  for  all  other  GPS  satellites  on  these  days  have  similar  biases. 
The  bias  between  the  broadcast  clocks  and  the  JPL  high  rate  clocks  is  due  to  a 
bias  between  the  reference  ground  clock  used  in  the  estimation  of  the  JPL  high 
rate  clocks  and  GPS  system  time  [48].  Fortunately,  the  bias  is  absorbed  in  the 
estimate  of  the  receiver  clock  bias,  since  it  is  nearly  constant  for  all  satellites. 


3-5  Feb  1997  -  prn02  bias  =  72.33  m 


Figure  2.16  -  Broadcast  vs.  High  Rate  Clock  Biases  for  PRN  2 


3-5  Feb  1997  -  prn09  bias  =  71.91  m 


Figure  2.17  -  Broadcast  vs.  High  Rate  Clock  Biases  for  PRN  9 

2.2.9.2  Application  of  JPL  High  Rate  GPS  Clock  Estimates 


While  application  of  the  JPL  high  rate  clock  estimates  does  not  completely 
remove  the  effects  of  SA,  due  to  the  reasons  mentioned  above,  it  does  provide  a 
pessimistic  indication  of  the  results  attainable  with  SA  off.  When  processing  the 
4  Feb  1997  GPS/MET  data,  gaps  in  the  high  rate  clock  data  are  interpolated 
through  and  where  data  is  not  available  for  a  GPS  satellite,  the  satellite’s 
measurement  is  not  used.  Table  2.12  shows  results  when  the  high  rate  GPS  clock 
estimates  are  used  instead  of  those  broadcast  for  single  and  dual  frequency  and  for 
precise  and  broadcast  ephemerides.  As  expected,  the  case  where  the  high  rate 


GPS  clocks  are  used  with  dual  frequency  ionospheric  correction  and  precise  GPS 
ephemerides  provide  the  most  accurate  and  most  precise  results. 


Table  2.12  -  GEODE  GPS/MET  SA  Off  Results  Comparison 


GPS/MET  -  4  Feb  1997 

Error  Mean  (m) 

RMS  Error  (m) 

Single  Frequency 

R 

•  T ’  • 

c 

R 

I 

C 

3D 

Broadcast  Clock  Biases 

Broadcast  GPS  Ephemerides 

0.55 

-0.99 

-0.10 

2.74 

10.24 

2.15 

10.81 

B  roadcast  Clock  B  iases 

Precise  GPS  Ephemerides 

0.58 

-1.05 

-0.09 

2.99 

10.11 

1.97 

10.73 

High  Rate  Clock  Biases 

Broadcast  GPS  Ephemerides 

0.28 

-3.59 

-0.04 

2.22 

7.23 

0.77 

7.60 

High  Rate  Clock  Biases 

Precise  GPS  Ephemerides 

0.30 

-3.56 

-0.03 

2.31 

6.95 

0.74 

7.36 

•  ■  ' "-V irirr.'r.V '.v 

Dual  Frequency 

R 

I 

C 

R 

I 

c 

3D 

Broadcast  Clock  Biases 

Broadcast  GPS  Ephemerides 

0.11 

2.04 

-0.11 

2.48 

9.68 

2.12 

10.21 

Broadcast  Clock  Biases 

Precise  GPS  Ephemerides 

0.11 

2.01 

-0.11 

2.48 

9.22 

2.03 

9.76 

High  Rate  Clock  Biases 

Broadcast  GPS  Ephemerides 

-0.08 

-1.38 

-0.04 

2.10 

6.14 

0.86 

6.55 

High  Rate  Clock  Biases 

Precise  GPS  Ephemerides 

-0.06 

-1.35 

-0.03 

1.91 

5.43 

1.04 

5.85 

2.3  SA  Free  T/P  Data 

Another  data  set  used  in  the  research  performed  for  this  dissertation  is 
from  the  T/P  satellite  (collected  on  5  May  2000).  Again,  SA  was  turned  off  on  2 
May  2000.  Additional  information  regarding  T/P  can  be  found  at: 
http://topex-www.jpl.nasa.gov/ 

While  the  dynamic  environment  and  ionospheric  effects  on  T/P  (1336  km 
altitude)  are  more  benign  than  those  at  lower  altitudes,  SA  free  T/P  data  is  also 
used  to  gauge  GEODE’s  performance. 


Figure  2.18  shows  plots  of  the  radial,  in-track  and  cross-track  errors 
between  GEODE’s  position  estimates  for  T/P  and  JPL’s  POE.  The  1.2  m  3D  RSS 
error  is  a  significant  improvement  over  the  GPS/MET  results  (6.55  m  3D  RSS) 
with  high  rate  GPS  clock  estimates  used  to  reduce  the  effects  of  SA.  Another 
indication  of  the  significance  of  SA  being  turned  off  is  the  measurement  residual 
plot  in  Figure  2.19.  A  comparison  between  a  plot  of  the  measurement  residuals 
for  the  GPS/MET  data  with  high  rate  GPS  clock  estimates  applied  (shown  in 
Figure  2.20)  and  the  T/P  measurement  residuals  (shown  in  Figure  2.19)  shows  an 
order  of  magnitude  difference.  The  T/P  measurement  residual  RMS  is  3.3  m 
while  the  GPS/MET  measurement  residual  RMS  is  31.5  m.  Table  2.13  shows  a 
summary  of  T/P  results  with  and  without  precise  GPS  ephemerides.  Improvement 
to  GEODE  will  be  gauged  in  this  dissertation  by  processing  both  GPS/MET  and 
T/P  data. 
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Figure  2.18  -  GEODE  SA  Free  T/P  Results 


Time  (hours) 


Figure  2.20  -  GPS/MET  with  High  Rate  GPS  Clock  Estimates  -  Measurement 

Residuals 


Table  2.13  -  SA  Free  T/P  Results  with  GEODE 


5  May  2000 

Error  Mean  (m) 

Error  RMS  (m) 

R 

i 

C 

R 

I 

C 

3D 

T/P 

0.21 

-0.04 

0.01 

0.35 

1.06 

0.44 

1.20 

TOEX  with  Precise  GPS 
Ephemerides 

0.22 

-0.16 

0.00 

0.34 

0.97 

0.31 

1.08 

2.4  Summary 

The  data  collected  by  both  the  GPS/MET  experiment  on  board  the 
OrbView-1  satellite  and  T/P  is  extremely  valuable  in  gauging  the 
accuracy/precision  of  the  real-time  satellite  OD  software  suite  known  as  GEODE. 
Several  GEODE  processing  schemes  are  investigated  and  results  are  presented  for 
several  modifications  to  GEODE  that  improve  performance  and  provide  insight 
into  GEODE’s  abilities  in  an  SA  free  GPS  environment.  Results  presented  in 
Table  2.12  and  Table  2.13  will  be  used  as  a  benchmark  to  compare  results 
achieved  with  suggested  improvements  to  GEODE  in  Chapters  3-6. 


CHAPTER  3 


Earth  Gravity 

Post-processing  and  near  real-time  Precise  OD  (POD)  systems  have  the 
luxury  of  using  state-of-the-art  work  stations  or  networked  systems,  while  real¬ 
time  systems  are  limited  to  the  latest  space  qualified  hardware.  Unfortunately,  the 
processing  capability  difference  between  ground  computers  and  those  found  on 
satellites  is  significant.  High-end  flight  computers  have  roughly  the  equivalent 
computational  power  of  a  60  MHz  Pentium  I  [49].  Therefore,  on  board  OD 
systems  must  minimize  the  computational  burden  of  their  software  while 
maximizing  accuracy  and  autonomy.  One  area  where  computational  burden  can 
be  reduced  significantly,  without  overly  decreasing  accuracy,  is  the  degree  and 
order  of  the  gravitational  acceleration  model  used. 

This  chapter  presents  the  results  of  an  investigation  into  the  effects  on 
accuracy  and  computational  burden  of  the  degree  and  order  of  the  gravity  model 
used  in  propagating  an  orbit  and  in  performing  satellite  OD.  Since  20%  of 
GEODE’s  computational  time  is  spent  evaluating  gravitational  accelerations,  a 
method  is  also  investigated  to  significantly  reduce  GEODE’s  computational 
burden  by  replacing  spherical  harmonic  coefficient  with  an  approximation 


method. 
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3.1  Model  Size  Integration  Results  and  Computational  Burden  Estimates 

The  size  (degree  and  order)  of  the  gravity  model  used  in  propagating  a 
satellite  orbit  from  one  time  to  the  next  significantly  affects  the  accuracy  of  the 
propagation.  But,  as  the  degree  and  order  of  the  gravity  model  increase,  so  does 
the  computational  burden  on  the  computer  used  for  the  propagation.  In 
performing  real-time,  autonomous  satellite  OD,  it  is  highly  desirable  to  minimize 
the  computational  burden  while  maximizing  accuracy.  The  accuracy  required  by 
the  integrator  depends  on  the  altitude  of  the  satellite  and  the  type  of  estimation 
filter  used.  The  higher  the  satellite,  the  smaller  the  size  of  the  gravity  model 
needed  to  achieve  the  same  accuracy.  Figure  3.1  shows  a  comparison  of 
accelerations  calculated  at  various  altitudes  between  a  full  70x70  JGM-2  gravity 
model  and  truncated  models.  The  trends,  as  expected,  show  decreasing  accuracy 
with  smaller  models  and  lower  altitude. 

The  filter  type  also  affects  the  size  of  the  gravity  model  needed.  If  the 
filter  is  a  batch  or  batch-sequential,  then  the  state  is  not  measurement  updated 
after  each  measurement  is  processed.  Therefore,  the  integration  time  used  in  a 
batch  or  batch-sequential  processor  is  typically  a  day  or  longer.  On  the  other 
hand,  an  EKF  updates  the  state  every  time  a  measurement  is  processed,  thereby 
significantly  reducing  the  integration  time.  In  the  case  of  processing  the 
GPS/MET  data,  the  measurements  are  10  seconds  apart  so  the  integration  time 
between  measurements  is  only  10  seconds. 
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Figure  3.1  -  Acceleration  Accuracy  Comparison  of  Various  Size  Gravity  Models 

at  Various  Altitudes 


To  accurately  assess  the  size  of  gravity  model  to  be  used  in  an  on  board 
EKF,  the  filter  must  be  run  with  real  or  simulated  data  using  various  size  gravity 
models.  The  position  and  velocity  estimates  are  then  compared  to  “truth”  to 
assess  accuracy. 

Another  method  that  can  be  used  to  guide  the  choice  of  the  size  of  the 
gravity  model  used  in  real-time  OD  is  to  perform  a  one-day  integration  using  an 
orbit  similar  to  the  satellite’s  mission  orbit.  This  process  can  be  used  to  narrow 
the  design  space  used  for  optimizing  accuracy/precision  versus  computational 
burden.  Table  3.1  shows  a  comparison  between  GPS/MET  orbit  propagations 


performed  with  several  truncations  of  the  EGM-96  gravity  model.  The  first 
column  is  the  number  of  coefficients  in  the  spherical  harmonic  coefficient 
expansion  used  to  calculate  the  acceleration.  The  second  column  shows  the 
number  of  Floating  Point  Operations  (FLOPS)  required  in  Matlab  to  calculate  the 
non-spherical  gravitational  acceleration  at  one  point.  The  third  column  shows  the 
time  required  to  propagate  a  low  Earth,  circular  orbit  for  one  full  day  using 
GEODE.  Each  of  these  estimates  of  computational  burden  are  roughly 
proportional,  as  expected.  The  final  column  shows  the  3D  RSS  error  of  the 
truncated  models  compared  to  the  full  70x70. 

The  initial  conditions  for  the  propagation  using  the  70x70  gravity  model 
are  shown  in  Table  3.1.  The  initial  conditions  of  the  propagations  performed  with 
truncated  gravity  models  are  adjusted  to  provide  the  best  fit  to  the  70x70  model. 


Table  3.1  -  Computational  Burden  and  Accuracy  of  Various  Truncations  of  the 

EGM-96  Gravity  Model 


Model 

Size 

Number  of 
Coefficients 

Number  of 
Matlab  Flops 

GEODE 
Integration 
Time  (sec) 

EGM-96 
RSS  Error 
Gompared  to 
70x70 (m) 

5x5 

36 

528 

51.5 

91.80 

10x10 

126 

1553 

52.1 

31.40 

20x20 

456 

5103 

55.8 

6.13 

30x30 

986 

10653 

64.8 

1.61 

40x40 

1716 

18211 

68.1 

0.61 

50x50 

2646 

27761 

77.9 

0.13 

70x70 

5106 

52853 

98.0 

N/A 

Table  3.2  -  Initial  Conditions  of  Propagation  With  EGM-96  70x70  Gravity  Model 


Osculating  Elements 

a 

7122.463  km 

e 

0.00143 

i 

70.00° 

Q 

196.54° 

G) 

233.66° 

o 

2.47° 

The  purpose  of  Table  3.1  is  to  show  the  computational  burden  required 
and  propagation  accuracy  attainable  with  truncations  of  the  EGM-96  model.  As 
shown  in  Table  3.1,  there  is  only  a  4.5  m  difference  in  accuracy  between  the 
20x20  truncation  and  the  30x30  truncation  but  over  a  16%  computational  burden 
difference.  With  knowledge  of  the  accuracy  attainable  with  various  truncations 
the  next  step  is  to  compare  OD  accuracy. 


3.2  GEODE  Results  and  Computational  Burden  Changes  with  Various 
Gravity  Model  Truncations 

There  is  a  moderate  improvement  in  the  24-hour  propagation  accuracy 
between  a  20x20  truncation  and  a  30x30  truncation  (roughly  4.5  m)  at  the 
GPS/MET  altitude  and  inclination.  This  section  presents  a  comparison  of  the 
GEODE  OD  accuracy  and  computational  burden  using  various  truncations  of  the 
JGM-2  and  EGM-96  models.  Table  3.3  presents  the  JGM-2  Earth  gravity  state 
noise  correlation  times  generated  with  Fortran  77  code  produced  by  GSFC  called 
autocor4.for  (see  Section  2.2.2).  Table  3.4  shows  GEODE’s  results  with  various 
truncations  of  the  JGM-2  gravity  model  using  the  data  in  Table  3.3. 


There  is  only  a  one-meter  difference  in  OD  accuracy  between  the  20x20 
truncation  and  the  higher  order  truncations.  The  best  results  are  obtained  using 
the  26x26  truncation.  It  is  also  interesting  to  note  that  although  there  is  only  a 
one-meter  accuracy  improvement  from  the  20x20  truncation  to  the  30x30 
truncation,  there  is  an  11%  difference  in  the  processing  time.  In  general,  the 
accuracy  and  precision  increase  as  the  degree  and  order  of  the  gravity  model 
increase. 

The  EGM-96  Earth  gravity  state  noise  correlation  times  are  listed  in  Table 
3.5  and  the  GEODE  results  are  listed  in  Table  3.6.  GEODE  performed  similarly 
with  both  gravity  models.  A  comparison  between  the  3D  position  error  RSS  for 
both  the  JGM-2  and  EGM-96  models  can  be  seen  graphically  in  Figure  3.2.  The 
conclusion  from  these  results  is  that  a  30x30  gravity  model  provides  the  best 
balance  of  accuracy,  precision  and  computational  burden.  Next  a  method  will  be 
presented  to  replace  GEODE’s  current  method  of  calculating  Earth  gravity  in  an 
attempt  to  maintain  or  improve  accuracy  and  precision  while  reducing  the 
computational  burden. 


Table  3.3  -  Earth  Gravity  JGM-2  State  Noise  Correlation  Times 


autocor4.for 

Direction 

R 

i 

c 

5x5 

386.2048 

24.5966 

754.0956 

10x10 

211.0938 

4.5542 

420.7908 

20x20 

136.1409 

0.1918 

271.4699 

21x21 

132.9159 

0.2454 

264.8726 

22x22 

129.8031 

0.1773 

258.6294 

23x23 

127.0829 

0.2589 

253.0201 

24x24 

124.4671 

0.2155 

247.7713 

25x25 

122.0885 

0.1718 

242.9791 

26x26 

119.9174 

0.1766 

238.5807 

:  27x27 

117.9191 

0.1658 

234.5325 

28x28 

116.1295 

0.1952 

230.8798 

29x29 

114.4664 

0.1525 

227.5483 

30x30 

112.9898 

0.15918 

224.5475 

40x40 

104.9078 

0.1340 

208.1667 

50x50 

102.8627 

0.1311 

204.0012 

70x70 

102.4124 

0.1301 

203.0799 

Table  3.4  -  GEODE  GPS/MET  JGM-2  Comparison 


GPS/MET 

4  Feb  1997 

..  Mean  Error  (m) 

RMS  Error  (m) 

Gravity 

Model 

RunTime 

(sec) 

R 

I 

C 

R 

'■  'I::. 

C 

3D 

05x05 

60.0 

1.71 

-5.57 

-0.06 

25.42 

53.87 

29.42 

66.43 

10x10 

64.0 

0.54 

-2.03 

0.62 

7.66 

18.16 

7.15 

20.96 

20x20 

64.4 

0.52 

-0.97 

0.02 

3.12 

11.05 

2.93 

11.85 

22x22 

68.7 

0.66 

-1.14 

0.05 

2.90 

10.53 

2.26 

11.15 

24x24 

69.0 

0.58 

-0.91 

-0.07 

2.77 

10.56 

2.56 

11.21 

25x25 

69.5 

0.58 

-0.88 

-0.07 

2.78 

10.36 

2.03 

10.92 

26x26 

71.1 

0.60 

-1.02 

-0.06 

2.69 

10.18 

2.06 

10.73 

27x27 

71.7 

0.57 

-0.97 

-0.06 

2.72 

10.26 

2.13 

10.82 

28x28 

72.0 

0.56 

-0.99 

-0.06 

2.74 

10.31 

2.26 

10.91 

30x30 

72.5 

0.55 

-0.99 

-0.10 

2.74 

10.24 

2.15 

10.81 

40x40 

76.9 

0.52 

-0.90 

-0.07 

2.71 

10.14 

2.29 

10.75 

50x50 

86.0 

0.52 

-0.88 

-0.05 

2.71 

10.12 

2.41 

10.75 

70x70 

117.8 

0.52 

-0.88 

-0.05 

2.70 

10.11 

2.38 

10.74 

Table  3.5  -  Earth  Gravity  EGM-96  State  Noise  Correlation  Times 


autocor4.for 

Direction 

R 

I 

C 

5x5 

387.2046 

25.1933 

754.8398 

10x10 

211.0605 

4.8500 

420.5247 

20x20 

134.3853 

0.3476 

267.6209 

21x21 

131.0949 

0.3915 

261.1110 

22x22 

128.4610 

0.3233 

255.8411 

23x23 

126.1396 

0.3875 

251.0632 

24x24 

123.9387 

0.3443 

246.6534 

25x25 

121.9078 

0.3012 

242.5654 

26x26 

120.0385 

0.3004 

238.7804 

27x27 

118.2884 

0.2866 

235.2359 

28x28 

116.7078 

0.3085 

232.0099 

29x29 

115.2452 

0.2680 

229.0789 

30x30 

113.9493 

0.2710 

226.4442 

40x40 

106.6019 

0.2344 

211.5380 

50x50 

104.5926 

0.2283 

207.4407 

70x70 

104.1469 

0.2267 

206.5280 

Table  3.6  -  GEODE  GPS/MET  EGM-96  Comparison 


GPS/MET 

4  Feb  1997 

Error  Mean  (m) 

RMS  Error  (m) 

Gravity  Model 

R 

.  I  . 

C 

R 

I  ; 

C 

3D 

i  05x05 

1.71 

-5.58 

-0.06 

25.45 

53.88 

29.40 

66.45 

10x10 

0.54 

-2.01 

0.62 

7.70 

18.21 

7.16 

21.03 

20x20 

0.52 

-0.96 

0.01 

3.13 

11.03 

2.94 

11.83 

22x22 

0.57 

-0.87 

-0.08 

2.78 

10.55 

2.63 

11.22 

24x24 

0.57 

-0.87 

-0.08 

2.78 

10.55 

2.63 

11.22 

25x25 

0.57 

-0.82 

-0.07 

2.82 

10.40 

2.11 

10.98 

26x26 

0.60 

-0.98 

-0.07 

2.71 

10.21 

2.12 

10.78 

27x27 

0.57 

-0.95 

-0.07 

2.73 

10.26 

2.16 

10.84 

28x28 

0.56 

-0.95 

-0.07 

2.75 

10.33 

2.26 

10.92 

30x30 

0.55 

-0.97 

-0.11 

2.75 

10.25 

2.18 

10.83 

40x40 

0.53 

-0.90 

-0.07 

2.71 

10.16 

2.35 

10.78 

50x50 

0.52 

-0.88 

-0.06 

2.71 

10.13 

2.42 

10.76 

70x70 

0.53 

-0.88 

-0.06 

2.71 

10.12 

2.40 

10.75 
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Figure  3.2  -  JGM-2  /  EGM-96  3D  Position  Error  RSS  Comparison 


3.3  Gravitational  Acceleration  Approximation  Function  (GAAF) 

Improving  the  speed  of  gravitational  acceleration  computations  can 
significantly  reduce  the  computational  burden  of  an  on  board  computer 
performing  real-time  satellite  OD.  Several  methods  currently  used  to  improve  the 
speed  of  computing  gravitational  accelerations  are  truncating  the  gravity  model, 
pre-selecting  a  significant  subset  of  coefficients,  or  tuning  a  truncated  gravity 
field  through  estimation.  Another  method  described  in  Hujsak  [50]  represents 
gravitational  accelerations  in  terms  of  an  instantaneous  two-body  acceleration  for 
an  instantaneous  pseudocenter  of  the  Earth.  A  set  of  pseudocenters  at  various 
heights  and  fixed  latitude  and  longitude  are  fit  to  polynomials  or  polynomial 
quotients.  The  polynomial  coefficients  are  calculated  at  various  latitudes  and 
longitudes.  The  stored  coefficients  are  used  to  calculate  the  gravitational 


66 


acceleration.  First,  polynomials  are  evaluated  using  the  satellite’s  height  to 
recover  a  set  of  six  pseudocenters  surrounding  the  desired  latitude  and  longitude 
of  the  satellite.  Then,  a  six-point,  bi-variate  interpolation  scheme  is  implemented 
to  recover  the  pseudocenter  at  the  given  latitude  and  longitude.  The  pseudocenter 
is  then  used  to  recover  the  gravitational  acceleration. 

In  the  following  development  of  GAAF,  the  method  of  pseudocenters  is 
introduced  and  the  benefits  of  ignoring  the  C2,o  geopotential  coefficient  term  in 
pseudocenter  formulation  are  discussed.  Next  polynomial  representation  of 
pseudocenters  as  a  function  of  height  is  shown,  followed  by  interpolation  of 
pseudocenters  on  a  sphere  of  common  height.  Finally,  storage  requirements  are 
shown. 

3.3.1  Method  of  Pseudocenters 

Given  the  Earth  Centered,  Earth  Fixed  (ECEF)  gravitational  acceleration 
on  a  spacecraft,  aECEF,  and  using  the  restricted  two-body  equation  of  motion,  a 
pseudocenter,  cECEF,  can  be  calculated.  In  equation  (3.1)  pECEF  is  a  pseudoradius 
to  the  spacecraft. 

aECEF=-F^fL  C3-1) 

P 

In  the  above  equation  aECEF  is  the  nonspherical  acceleration  calculated 
using  spherical  harmonic  coefficient  expansion  of  the  nonspherical  gravitational 
coefficients.  Now  the  actual  radius  to  the  satellite  is  the  pseudoradius  plus  a 
pseudocenter.  A  pseudocenter  is  a  vector  from  the  center  of  the  Earth  to  where 
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the  center  of  the  Earth  would  need  to  be  for  the  nonspherical  acceleration  acting 
on  the  satellite  to  equal  the  restricted  two-body  acceleration  acting  on  the  satellite. 

*ECEF  =  PeCEF  +  CECEF  (3-2) 


or 


CECEF  —  %CEF  PeCEF 


By  definition. 


and 


P  =  |p 


ECEF 


(3.3) 


(3.4) 


Pecef-  P  Pecef  (3-5) 

Also,  since  the  pseudoradius  vector  is  in  the  opposite  direction  of  the  acceleration: 


fr  _  aECEF 

Pecef  i— 


(3.6) 


•‘ECEF 


Therefore,  substituting  equation  (3.5)  into  equation  (3.3)  yields: 


CECEF  —  rECEF  P  PeCEF 


Now,  dot  product  both  sides  of  equation  (3.1)  with  p  yielding: 


i-  i_ _ L 

aECEF  —  P  2 

P 


Therefore, 


(3.7) 


(3.8) 


P=  fe 


P 


(3-9) 


‘ECEF 


Leading  to  the  result: 


-  P  ^ 

CECEF  —  rECEF  -  -  PeCEF 
V  aECEF 


(3.10) 
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or 


CECEF  “  rECEF  + 


ECEF 


|3/2 


ECEF  | 


(3.11) 


3.3.2  Benefits  of  Pseudocenter  Formulation  Without  the  C2,0  (J2)  Term 


As  noted  in  HujaK  [50]  the  magnitude  of  the  pseudocenters  are  less  than 
15  km  (|cECEF|  <  15km)  for  all  heights  above  100  km.  However,  omitting  the  C2,o 

term  from  the  gravitational  acceleration  calculation  yields  |cECEF|  <  250  meters 

which  can  halve  the  storage  requirements  for  the  support  data.  The  omitted  C2,o 
acceleration  is  then  added  in  after  the  pseudocenter  is  used  to  recover  the 
acceleration.  Let 


aECEF0  —  aECEF  aECEFC2  0 


(3.12) 


where 


ECEF<- 


C2,0 


r  lEL 
2’°  2  r3 


Rr 


\2 


V  r  J 


5|  — 


-1 


3  py  f 
2,0  2  r3  I  r 


51  — 


V  J 

2(  /  \2  \ 

-1 

r) 


v 


3  pzf  R 


_C 

2,0  2  r3 


_E 

V  r  J 


\2(  fz'2  ^ 


51 7  -3 


V 


[51] 


and  C2l0  =  -h  =  -1.082626925638815xl0'03  (from  JGM-2) 
Now  the  pseudocenter  is  calculated  independent  of  C2,o* 


“'ECEFo 


=  r, 


ECEF 


AECEF0 


3/2 

aECEF0 

(3.13) 


(3.14) 
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3.3.3  Polynomial  Representation  of  Pseudocenter  as  a  Function  of  Height 


The  first  step  in  creating  a  set  of  coefficients  by  which  gravitational 
acceleration  can  be  recovered  is  to  determine  the  altitude  range  of  the  spacecraft. 
Hujsak  [50]  uses  an  altitude  range  of  400  km  to  1500  km.  However,  since  the 
proposed  real-time  OD  system  will  be  used  with  spacecraft  in  circular  or  near 
circular  orbits,  smaller  ranges  are  investigated  here.  Hujsak  [50]  suggests  using 


quotients  of  polynomials  of  the  form: 

_  a0  +  atx  +  a2x2  +. . .  .+an_1xn_1 

1  +  bjX+. . .  .+bd_jXd 


(3.15) 


where 


h-hmm 

U  _ 1_ 

11  max  11  min 


(3.16) 


Hujsak  [50]  reports  n  +  d  <  8  if  the  height  in  equation  (3.16)  is  limited  to 
hmax  =  1500  km  >  h  >  h„un  =  400  km. 

In  this  study,  the  GPS/MET  spacecraft,  altitude 
hmax  =  752  km  >  h  >  h^n  =  732  km, 

is  considered  first.  Up  to  seventh  order  polynomials  of  the  form 

C;  =  a0 +ajX  +  a2x2+...+a7x7  (3.17) 

are  considered  for  the  ease  of  implementation. 
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3.3.4  Pseudocenter  Interpolation  on  a  Sphere  of  Common  Height 

Hujsak  [50]  suggests  a  six-point,  bi-variate  interpolation  scheme  to  take 
the  pseudocenters  calculated  from  the  coefficients  to  the  pseudocenter  for  a  given 
latitude  and  longitude.  Figure  3.3  shows  the  points  used  in  the  six-point,  bi¬ 
variate  interpolation  scheme  [52],  The  circles  represent  the  pseudocenters 
calculated  at  the  fixed  latitudes  and  longitudes  while  the  square  represents  the 
latitude  and  longitude  of  the  current  satellite  position. 


Figure  3.3  -  Six-Point  Bi-Variate  Interpolation  Scheme 


The  equation  used  for  the  six-point,  bi-variate  interpolation  is: 
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CofopAq.h) =  ^o(4>0  +pAcf),A,0  +  qAA,,h) 
=.5q(q  - 1)  c(cj)0 ,  X0  -  AX,  h) 
+.5p(p-l)c(<|>0-A<|>,X0,h) 

+  (l  +  pq-p2-q2)c(<j)0A0,h) 
+.5p(p  -  2q  + 1)  c(<J)0  +  A<>,  X0 ,  h) 
+.5q(q  -  2p  + 1)  c(<|)0 ,  X0  +AX,h) 
+  pq  c(<)>0  +  A<|),  X0  +  AX,  h) 

where 

(<I>p — 4>o) 

p"  A4> 


(3.18) 


(3.19) 


3.3.5  GAAF  Storage  Requirements 

Hujsak  [50]  suggests  larger  increments  of  longitude  as  latitude  increases  to 
reduce  storage  requirements.  In  order  to  simplify  indexing  into  the  arrays  holding 
the  pseudocenter  coefficients,  there  is  an  overlap  between  sets  of  pseudocenter 
coefficients,  i.e.,  latitudes  -59°  and  -60°  are  contained  in  sets  3b  and  2b  (see  Table 
3.7).  The  longitudes  contained  in  3b  are  -2°,  0°,  2°,.. .,360°.  Single  precision  (4 
byte)  floating  point  variables  are  used.  The  number  of  storage  bytes  needed  for 
each  latitude  and  longitude  is  calculated  by: 

#  coefficients/pseudocenter  *  3  pseudocenters  *  4  bytes 


At  the  GPS/MET  altitude,  1,774,296  bytes  are  required  (1.69  MB) 
compared  to  5,024,832  bytes  (4.79  MB)  for  the  original  implementation.  The 
tradeoff  is  a  smaller  altitude  and  inclination  range  but  the  storage  reduction  is 
significant.  Table  3.7  shows  the  longitude  increments  for  different  latitudes  and 


the  storage  requirements  for  each  range. 

Table  3.7  -  Original  GPS/MET  GAAF  Storage  Requirements  For  3  Coefficients 


latitudes 

AX. 

#lats 

#lons 

#bytes  each 

total  bytes 

3b 

-75°  to  -59° 

2.0° 

17 

182 

36 

111,384 

2b 

-61°  to -47° 

1.5° 

15 

242 

36 

130,680 

1 

-49°  to  49° 

1.0° 

99 

362 

36 

1,290,168 

2a 

47°  to  61° 

1.5° 

15 

242 

36 

130,680 

3a 

59°  to  75° 

2.0° 

17 

182 

36 

111,384 

Total 

1,774,296 

To  decrease  the  complexity  of  implementing  GAAF,  coefficients  are  also 
generated  using  a  single  longitude  increment  of  1°.  The  number  of  bytes  required 
to  store  the  simplified  implementation  is  151  *  362  *  3  *  3  *  4  or  1,967,832  bytes 
or  (1.88  MB).  The  small  increase  in  the  number  of  bytes  required  is  worth  the 
improved  simplicity  of  implementation.  A  detailed  explanation  of  how  GAAF  is 
implemented  and  validated  is  shown  next. 


3.3.6  Integration  With  GAAF 

The  first  task  in  implementing  GAAF  is  to  create  a  capability  to  generate 
highly  accurate  gravitational  accelerations  at  multiple  latitudes,  longitudes  and 
heights.  Since  GEODE  already  has  a  JGM-2  30x30  gravity  model  working, 
“subroutines”  (classes)  are  taken  from  GEODE  and  augmented  to  include  the  full 
JGM-2  70x70  gravity  model.  The  updated  subroutines  are  then  reintroduced  into 


GEODE.  This  is  done  to  validate  the  new  gravity  model  and  to  provide  a  method 
for  evaluating  GAAF.  A  GPS/MET  orbit  is  integrated  with  GEODE  and  the 
results  are  compared  to  an  integration  performed  with  Analytical  Graphics’ 
Satellite  Toolkit  (STK).  Additional  information  about  STK  can  be  found  at 
http://www.stk.com.  A  direct  comparison  in  the  Earth  Centered  Inertial  (ECI) 
coordinate  system  shows  a  3D  RSS  error  of  0.738  meters.  There  is  no  definitive 
explanation  for  the  difference,  except  possibly  differences  in  the  way  both 
programs  handle  polar  motion,  differences  in  nutation  and  precession  constants 
and  differences  in  the  UTC  to  UT1  time  conversion.  Figure  3.4  shows  the 
difference  between  the  70x70  GEODE  trajectory  and  the  70x70  STK  trajectory. 

3D  RSS  Error  =  0.738  m 


Figure  3.4  -  GEODE  /  STK  Integration  Comparison 
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Another  method  used  to  ensure  the  new  70x70  gravity  model  is  working 
properly  is  to  compare  the  difference  between  a  70x70  and  30x30  propagation  in 
GEODE  and  the  difference  between  a  70x70  and  30x30  propagation  in  STK.  The 
“differences”  comparison  shows  that  truncating  the  gravity  model  manifests  itself 
the  same  way  for  both  the  STK  and  GEODE  propagations.  Figure  3.5  shows  the 
excellent  agreement  between  the  STK/GEODE  truncation  comparisons. 


The  conclusion  from  these  results  is  the  JGM-2  70x70  gravity  model 
implemented  in  GEODE  is  working  correctly.  The  JGM-3  and  EGM-96  gravity 
models  are  implemented  and  validated  in  GEODE  in  the  same  way. 
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The  next  step  is  to  implement  GAAF  with  the  GEODE  integrator  and 
compare  the  results  to  a  70x70  propagation.  With  the  gravitational  accelerations 
now  available,  a  method  is  devised  to  determine  the  height  range,  height 
increment  and  number  of  polynomial  coefficients  to  use  for  the  GPS/MET  orbit. 

Since  the  OrbView-1  satellite  carrying  GPS/MET  is  in  a  near  circular  orbit  with  a 
maximum  altitude  of  just  over  749  km  and  a  minimum  altitude  of  just  over  734 
km,  an  altitude  range  of  732  -  752  km  is  chosen.  The  decimated  range  is  chosen 
to  attempt  to  reduce  the  size  of  the  memory  needed  to  store  the  GAAF  data.  2  to 
8  polynomial  coefficients  are  used  to  fit  the  pseudocenters.  Hujsak  [50]  suggests 
using  a  quotient  of  polynomials  as  in  equation  (3.15)  but  this  method  could  not  be 
duplicated,  so,  polynomials  are  used  as  shown  in  equation  (3.17). 

Pseudocenters  at  fixed  latitude  and  longitude  are  fit  to  polynomials  using 
the  following  methodology: 

Example: 

let  height  increment  for  fit  Ah  =  2  km 
hmin  =  732  km,  lw  =  752  km 
<j)  =  30° 

X  =  30° 

n  s  number  of  pseudocenters  to  fit  =  1 1 
p  =  number  of  polynomial  coefficients  =  2  through  8 

therefore  from  equation  (3.16) 

x  =  [0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.0]T  (3.20) 

a  =  [a0  a,  a2  ...  ap_,] 


(3.21) 
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Now  the  normal  equation  can  be  defined  as: 

y  =  Xa+ ' 
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(3.22) 


(3.23) 


(3.24) 


where  the  solution  to  the  normal  equation  for  a  yields  the  coefficients  of  the 
polynomial  that  “fits”  the  vector  of  pseudocenters,  y .  There  are  many  ways  to 
solve  the  normal  equation.  The  method  of  least  squares,  as  detailed  in 
Montgomery  [53],  is 

a  =  (XTX)_1XTy  (3.25) 


However,  due  to  the  computational  burden  of  performing  an  inverse  on  a  matrix 
as  large  as  50x50  (a  much  larger  matrix  is  required  for  fitting  a  larger  height 
range),  a  Givens  orthogonal  transformation  is  used  to  solve  for  a  [39].  Here  an 
orthogonal  transformation  is  performed  to  upper  triangularize  the  matrix 

[X  y]  (3.26) 


resulting  in  a  matrix  of  the  form 
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0 

0 

0 


...  hln  b, 

0  hn"  > 

0  0 


(3.27) 


Then,  back  substitution  is  performed  to  solve  for  a ,  where  an  =  bn/hnn  and  so  on. 
Orthogonal  transformations  require  less  computational  burden  to  solve  the  normal 
equation  [39].  Table  3.8  shows  the  pseudocenter  elements  at  <(>  =  30°  and  A,  =  30° 
for  the  height  range  of  732-752  km  at  2  km  increments.  Figure  3.6  shows  the 
linear  nature  of  the  pseudocenters. 


Table  3.8  -  Test  Pseudocenters 


height 

X 

yi 

:  y3 

:  732 

0 

26.1723 

-40.4541 

32.0436 

734 

0.1 

26.2094 

-40.5061 

32.0783 

736 

0.2 

26.2465 

-40.5578 

32.1126 

738 

0.3 

26.2835 

-40.6093 

32.1466 

4.  740 

0.4 

26.3204 

-40.6605 

32.1803 

4  742 

0.5 

26.3572 

-40.7115 

32.2135 

744 

0.6 

26.3940 

-40.7623 

32.2465 

746 

0.7 

26.4307 

-40.8128 

32.2790 

748 

0.8 

26.4674 

-40.8631  j 

32.3113 

750 

0.9 

26.5040 

-40.9131 

32.3432 

752 

1 

26.5405 

-40.9629 

32.3747 
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Figure  3.6  -  Linear  Nature  of  Pseudocenters 


Next,  the  most  appropriate  number  of  polynomial  coefficients  to  fit  the 
pseudocenters  is  determined.  The  Residual  Mean  Square  (RMS)  (o2)  [53]  is 
calculated  for  set  of  polynomials  where 

2  yTy-aTXTy 

<7  = - 

n-p 


(3.28) 


Table  3.9  shows  the  polynomial  coefficients  calculated  and  their 
corresponding  residual  mean  square  error.  At  30°  latitude,  30°  longitude  and  with 


a  2  km  height  increment,  the  best  fit  is  p  =  3.  To  determine  the  best  height 
increment  and  number  of  polynomial  coefficients  for  the  GPS/MET  orbit,  a  range 
of  height  increments  from  0.25  to  4  km  are  used  along  with  a  range  of  polynomial 


coefficients  from  2  to  8  (first  order  to  seventh  order).  The  RMS  error  is 
calculated  at  multiple  (100’s)  latitudes  and  longitudes.  The  best  fit  is  determined 
by  the  average  of  the  RMS  errors.  The  GPS/MET  orbit  produced  a  best  fit  with 
the  height  increment  at  3.0  km  and  p  =  3.  An  example  set  of  coefficients  and 
RMS  values  for  various  order  polynomials  are  shown  below  in  Figure  3.7. 


Table  3.9  -  Example  GPS/MET  GAAF  Coefficients  and  RMS  for  yl 


a0 

al 

a2 

a3 

a4 

a5 

a6 

a7 

:■  cr2 

p  =  2 

26.173 

0.3681 

1.05e-07 

p  =  3 

26.172 

0.3714 

-0.00332 

2.80e-ll 

P  =  4 

26.172 

0.3714 

-0.003046 

-0.000182 

2.10e-ll 

p=5 

26.172 

0.3714 

-0.003036 

-0.00020 

8.36e-6 

9.35e-10 

p  =  6 

26.172 

0.3714 

-0.003037 

-0.000197 

5.93e-6 

9.71e-7 

1.45e-08 

p  =  7 

26.172 

0.3714 

-0.003038 

-0.000193 

-2.41e-6 

8.44e-6 

-2.50e-6 

5.06e-07 

P  =  8 

26.172 

0.3714 

-0.003040 

-0.000177 

-4.53e-5 

7.44e-5 

-4.95e-5 

1.25e-5 

7.08e-05 
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Figure  3.7  -  GPS/MET  Polynomial  Fit  RMS  Error 

The  pseudocenter  coefficients  are  then  calculated  and  stored  in  a  file 
formatted  as  an  array  in  the  C  programming  language.  The  array  variable  is  then 
initialized  as  a  static  variable  in  C  and  the  coefficients  are  accessed  according  to 
latitude  and  longitude.  Six  sets  of  coefficients  are  accessed  each  time  an 
acceleration  is  calculated.  Below  find  the  equations  used  to  recover  acceleration 
from  the  pseudocenter  coefficients  with  p  =  3.  The  six  pseudocenters  at  the  given 
height  are  calculated  to  interpolate  between  the  fixed  latitudes  and  longitudes. 
Equation  (3.18)  is  used  for  the  bi-variate  interpolation.  First,  the  pseudocenters 
are  calculated  from  the  tables  of  polynomial  coefficients. 


CECEF0  =  a0  +a,x  +  a2x2  i  =  1,2,3 
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(3.29) 


where  and  ao,  a\  and  &2  are  the  polynomial  coefficients. 


x  = 


h-h„ 


1_  _  U 

max  ~~  11  min 


(3.30) 


hmin  =  732  km 
hmax  =  752  km 

Next,  six  sets  of  pseudocenters  are  used  to  interpolate  using  equation  (3.18). 
Then  the  pseudoradius  is  calculated  by 


PeCEF  =  rECEF  —  CECEF0 

and  the  acceleration  (without  the  C2,o  term  is  calculated) 

—  _  ~M- Pecef 

“ECEF0  -  p3 

aECEF  =  aECEF0  +  aECEFC2  0 

where  aECEFc  is  found  by  equation  (3.13).  The  acceleration  in  the  ECEF 

coordinate  system  is  then  rotated  to  the  J2000  system  by  the  transpose  of  the 
transformation  matrix  [J2000toECEF].  Here 


(3.31) 


(3.32) 

(3.33) 


[  J2000toECEF]  =  [TODtoECEF]  [MODtoTOD]  [J2000toMOD]  (3.34) 
where  the  definition  of  the  transformation  matrices  above  can  be  found  in  Lee 


[37], 

With  the  code  written  to  recover  gravitational  accelerations  from  the  array 
of  polynomial  coefficients  a  24-hour  GPS/MET  propagation  is  performed  with  the 
GEODE  integrator.  However,  here  GAAF  is  used  instead  of  spherical  harmonic 
coefficient  expansion.  First  order  polynomials  are  used  to  fit  the  pseudocenters  to 


minimize  storage  requirements.  The  results  are  compared  to  a  propagation 
performed  using  the  30x30  and  70x70  JGM-2  models.  Table  3.10  shows  the  orbit 
errors  and  3D  RSS  error  of  the  difference  between  the  30x30  and  the  70x70  and 
GAAF  and  the  70x70.  The  run  times  for  the  various  integrations  are  also 
included  in  the  table.  The  GAAF  run  time  is  19.3%  shorter  than  the  30x30  model 
and  more  accurate.  GAAF  works  wonderfully,  increasing  accuracy  while 
decreasing  computational  burden.  Figure  3.8  shows  a  plot  of  the  propagation 
errors. 


Table  3.10  -  GPS/MET  GAAF  and  30x30  Integration  Compared  to  70x70 


GPS/MHT  -  4  Feb  1997 

”  :  :  Integration  RMS  Error  (m) 

Gravity  Model 

Run  Time  (sec) 

R 

I 

c 

3D 

70x70 

121.5 

30x30  vs  70x70 

74.4 

2.41 

1.03 

2.57 

3.67 

GAAF*  vs  70x70 

60.1 

0.14 

0.06 

0.13 

0.02 

*GAAF  created  using  70x70  JGM-2  Gravity  Model 
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Time  (hours) 


Figure  3.8  -  GPS/MET  GAAF  and  30x30  Integration  Compared  to  70x70 


3.4  GAAF  in  GEODE 

The  final  step  in  assessing  the  accuracy  and  computational  burden  of 
GAAF  is  to  use  it  in  GEODE  for  satellite  OD.  Table  3.11  shows  the  GEODE  OD 
results  using  GAAF.  The  results  in  Table  3.1 1  confirm  the  claim  in  Hujsak  [50] 
that  the  computational  burden  of  GAAF  is  equivalent  to  the  computational  burden 
of  a  5x5  gravity  model  with  the  accuracy  of  a  70x70.  GAAF  maintains  accuracy 
in  GEODE  and  significantly  reduces  the  computational  burden  with  first  order 
polynomials  (2  coefficients  per  pseudocenter  direction).  This  first  order 
GPS/MET  implementation  of  GAAF  requires  only  1.28  MB  of  storage  space. 


Table  3.11  -  GEODE  GPS/MET  Gravity  Model  Results  Comparison 


GPS/MET 

4  Feb  1997 

Mean  Error  (m) 

RMS  Error  (m) 

Gravity  Model 

*Run 

Time 

(sec) 

R 

I 

C 

R 

I 

c 

3D 

JGM-2  05x05 

39.5 

1.71 

-5.57 

-0.06 

25.42 

53.87 

29.42 

66.43 

JGM-2  30x30 

45.6 

0.55 

-0.99 

-0.10 

2.74 

10.24 

2.15 

10.81 

JGM-2  70x70 

83.2 

0.52 

-0.88 

-0.05 

2.70 

10.11 

2.38 

10.74 

GAAF 

39.2 

0.52 

-0.88 

-0.05 

2.71 

10.11 

2.38 

10.74 

EGM-96  30x30 

45.6 

0.55 

-0.97 

-0.11 

2.75 

10.25 

2.18 

10.83 

JGM-3  30x30 

45.6 

0.52 

-0.91 

-0.11 

2.74 

10.24 

1.98 

10.78 

Run  on  450MHZ  Pentium  n,  ] 

128  MB  RAM 

3.5  JGM-3  and  EGM-96  Gravity  Models 


The  JGM-3  and  EGM-96  gravity  models  are  also  implemented  in 
GEODE.  Table  3.11  shows  no  improvement  in  OD  accuracy  or  precision  in 


GEODE  when  different  models  are  used  instead  of  the  JGM-2  model. 


3.6  Summary 

Although  there  is  a  moderate  increase  in  the  propagation  accuracy 
between  a  JGM-2  20x20  gravity  model  and  a  JGM-2  30x30  gravity  model,  the 
OD  accuracy  and  the  computational  burden  difference  between  the  two  JGM-2 
models  is  minimal.  The  best  choice  in  balancing  accuracy/precision  and 
computational  burden  in  a  real-time  OD  system  appears  to  be  a  30x30  gravity 
model.  Also,  the  type  of  model  (JGM-2,  JGM-3  or  EGM-96)  does  not  appear  to 
affect  OD  accuracy  with  GEODE. 


A  different  approach  significantly  reducing  computational  burden  while 
maintaining  accuracy  is  Hujsak’s  Gravity  Acceleration  Approximation  Function 
(GAAF)  [50].  GAAF  maintains  the  accuracy  of  a  70x70  gravity  model  with  the 
computational  burden  of  a  5x5  model.  GAAF  improves  GEODE’s  computational 
burden  by  14%  compared  to  GEODE  with  a  30x30  gravity  model.  The  cost  of 
using  GAAF  is  a  1  to  5  MB  RAM  requirement,  depending  on  the  application. 
GAAF  should  definitely  be  considered  as  a  viable  alternative  to  conventional 
methods  of  calculating  the  gravitational  acceleration  in  real-time  OD  systems 
where  computational  burden  needs  to  be  minimized. 


CHAPTER  4 


Ionosphere 

Ionospheric  errors  can  significantly  degrade  the  accuracy  of  the  GPS 
pseudorange  measurement.  These  errors  can  range  from  a  few  meters  to  many 
tens  of  meters  and  with  SA  off,  ionospheric  errors  can  be  the  largest  measurement 
error  source  [10].  Dual  frequency  GPS  users  can  determine  ionospheric  errors,  to 
first  order,  using  a  combination  of  the  PI  and  P2  measurements  as  shown  in 
equation  (2.11).  Single  frequency  users  must  either  rely  on  ionospheric  models  or 
a  linear  combination  of  the  pseudorange  and  phase  measurements.  In  a  real-time 
OD  system  on  board  a  satellite,  the  overhead  of  using  a  model  may  be  unrealistic. 
Therefore,  the  Differenced  Range  Versus  Integrated  Doppler  (DRVID)  [29] 
method  is  presented  and  investigated  here.  The  results  of  a  DRVID 
implementation  in  GEODE  are  also  presented. 

4.1  Differenced  Range  Versus  Integrated  Doppler  (DRVID) 

The  ionosphere  is  the  region  of  the  upper  atmosphere  between  50  and 
1000  km  altitude.  It  contains  electrons  and  ions  formed  by  the  ionizing  radiation 
of  the  Sun.  These  electrons  and  ions  are  in  sufficient  quantities  to  significantly 
affect  radio  wave  propagation  where  the  delay  experienced  by  these  signals  is 
proportional  to  the  number  of  electrons  in  the  signal’s  path.  Unfortunately,  the 


affects  on  radio  wave  propagation  are  highly  variable  and  change  by  at  least  one 
order  of  magnitude  through  the  course  of  each  day.  The  effects  can  be: 

•  three  times  greater  on  a  signal  transmitted  near  the  horizon  compared  to 
one  transmitted  from  the  zenith  direction 

•  five  times  greater  during  the  day  than  at  night 

•  four  times  greater  in  November  than  July 

•  three  times  greater  during  solar  maximum  than  solar  minimum  [54] 

Because  of  this  high  variability,  the  ionosphere  is  very  difficult  to  model. 
There  are  four  regions  of  the  ionosphere  formed  by  different  chemical  interactions 
with  the  ultraviolet  (UV)  radiation  from  the  Sun.  Only  one  of  these  is  of  primary 
concern  for  satellite  OD  using  GPS.  It  is  the  F2  region  from  210  -  1000  km  and 
is  the  densest  and  has  the  highest  variability  [43].  The  peak  density  in  the  F2 
region  varies  between  250  and  400  km.  The  F2  region  can,  potentially,  cause  the 
most  significant  effects  on  GPS  receiving  systems  [43].  Another  area  associated 
with  the  ionosphere  is  the  protonosphere.  Its  region  designation  is  H+  and  it 
occupies  the  area  above  1000  km  and  extends  to  the  GPS  orbits.  It  is  composed 
of  ionized  hydrogen  and  some  helium  gas,  is  of  very  low  density  and  its  estimated 
contribution  to  signal  time  delay  is  10%  during  the  day  and  approximately  50%  at 
night  [43].  The  protonosphere  does  not  change  significantly  between  day  and 
night  but  is  depleted  during  major  magnetic  storms  and  can  take  several  days  to 
recover  [43]. 

The  ionosphere  can  be  a  significant  error  source  in  using  GPS 
measurements.  These  errors  can  range  from  a  few  meters  to  many  tens  of  meters 
at  the  GPS  frequencies.  There  are  effectively  eight  ways  the  ionosphere  effects 


GPS  measurements  [43],  but  only  two  are  significant  enough  to  be  discussed  here: 
group  delay  of  the  signal  modulation  (absolute  range  error)  and  carrier  phase 
advance  (relative  range  error). 

Most  receivers  (unless  codeless  or  Precise  Positioning  System  (PPS) 
capable)  can  only  collect  data  on  the  LI  frequency  and  thus  the  group  range  delay 
and  phase  advance  of  the  GPS  observations  cannot  be  calculated.  Modeling  the 
ionosphere  is  marginally  accurate  and  scaling  the  coefficients  to  satellite  altitudes 
may  cause  further  inaccuracies  [55].  Transmitting  measurements  or  delay 
corrections  to  a  satellite  introduces  additional  complexity  and  overhead. 
Unfortunately,  the  systematic  errors  introduced  by  not  modeling  the  ionospheric 
range  delay  and  phase  advance  can  hinder  the  estimation  of  unmodeled  or 
mismodeled  forces  using  techniques  like  Dynamic  Model  Compensation  (DMC) 
or  the  Reduced  Dynamic  Technique  (RDT)  [9].  In  Gold  [56]  it  is  shown  that 
RDT  performs  worse  than  a  purely  dynamic  run  if  ionospheric  errors  are  not 
removed.  Fortunately,  the  DRVID  measurement  can  potentially  remove  the 
systematic  ionospheric  delay  errors  that  can  reduce  the  accuracy  of  OD  systems 
estimating  empirical  accelerations. 

The  group  delay  is  proportional  to  the  Total  Electron  Content  (TEC)  and 
inversely  proportional  to  the  frequency  squared  of  the  modulated  signal.  The 
TEC  is  the  number  of  electrons  in  a  vertical  column  having  aim2  cross-section 
and  extending  from  the  receiver  to  the  GPS  satellite.  The  following  equation  is 
used  to  calculate  the  group  delay  in  units  of  time  [10]: 
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40  3 

Atp  =  — rTEC  (4.1) 

P|o°  cf 

where 

c  =  the  speed  of  light  in  a  vacuum  =  299,792,458  m/s 
f  =  signal  frequency  i.e.  fLi  =  1575.42  MHz  and  fL  =  1227.6  MHz 

or  expressed  in  units  of  range 

40  3 

APta,=^-TEc  (4.2) 

While  the  range  measurement  is  delayed,  the  ionosphere  advances  phase.  Phase 
advance,  in  units  of  time,  is  calculated  by 

40  3 

At^=-^TEC  (4.3) 

or  expressed  as  phase  range 

40  3 

AP^="— TEC  (4-4) 

Techniques  for  dealing  with  ionospheric  errors  include  modeling,  using 
direct  measurements,  using  a  dual  frequency  receiver  or  combining  range  and 
phase  measurements  on  a  single  frequency  receiver. 

There  are  many  models  available  to  estimate  TEC  to  attempt  to  calculate 
the  range  delay  and/or  phase  advance  of  GPS  measurements.  The  simplest  and 
most  readily  available  is  the  Klobuchar  model  consisting  of  eight  parameters 
contained  in  the  broadcast  navigation  message.  Surprisingly,  the  simple 
Klobuchar  model  consistently  exceeds  its  design  goal  of  50%  accuracy. 


Unfortunately,  very  complex  models  using  hundreds  of  coefficients  can  only 
marginally  improve  accuracy  and  cannot  consistently  surpass  75%  accuracy  [30]. 

Another  method  of  generating  ionospheric  delay  corrections  is  JPL’s 
Global  Ionosphere  Map  (GIM).  Dual  frequency  GPS  measurements  are  taken 
from  a  network  of  GPS  receivers.  These  measurements  are  processed  and  maps 
of  the  electron  content  of  the  ionosphere  are  produced.  GIM  will  be  integrated 
into  the  FAA’s  operational  WAAS  software  [16].  GIM  operations  can  take  place 
in  real-time,  on  board  a  satellite;  however,  there  is  significant  overhead  in 
providing  the  maps  to  the  satellite  for  processing.  It  is  unknown  what  the 
computational  burden  of  the  GIM  software  is  but  it  requires  a  383  KB  file  each 
day  [57].  WAAS  or  WADGPS  will  broadcast  ionospheric  corrections  but,  again, 
neither  system  is  currently  operational. 

The  DRVTD  technique  can  be  attributed  to  MacDoran  [29]  although 
previous  work  had  been  done  as  early  as  1966  [31].  Due  to  similarities  between 
the  original  signal  structures  and  those  used  in  GPS,  DRVID  is  generalized  for 
use  with  single  frequency  GPS  measurements  [58].  The  DRVID  technique  is 
further  investigated  in  Schreiner  [31]  and  Gold  [59].  JPL’s  implementation  of 
DRVID  in  satellite  OD  using  GPS  is  called  GRAPHIC  (Group  and  Phase 
Ionosphere  Calibration)  [30,  32,  60]. 

It  is  important  to  note  that  since  the  pseudorange  measurement  is  20-100 
times  less  precise  than  the  phase  measurement,  the  error  in  the  combined  DRVID 
measurement  is  half  that  of  the  pseudorange  measurement.  The  new 


91 

measurement  gets  the  worst  of  both  data  types  since  it  gets  the  reduced 
information  content  of  the  phase  delay  (due  to  the  phase  integer  ambiguity)  and 
the  approximate  noise  level  of  the  pseudorange  measurement  [60]. 


4.2  DR VID  Development 

The  geometric  range  from  the  user  satellite  to  each  GPS  satellite  is  given 


by: 


Pi  —  ^(Xu XGPS,  )  +(yu  yGPS;  )  +  (Zu  ZGPSj  )  (4.5) 

The  measured  code  based  pseudorange  in  meters  is: 

Ppi,  =  Pi  +  buc  -  bGPS.  c  +  Apion  (4.6) 

where: 


c  =  the  speed  of  light 

bu  =  user  clock  bias 

bops  =  GPS  spacecraft  clock  bias 

Apion  =  change  in  range  due  to  ionospheric  affects 

The  measured  beat  phase  based  pseudorange  in  cycles  is  [31,  54]: 

<t>Li,  +Nj  =-— Pj  -  buf  +  bGPS  f  +  -  Apion  (4.7) 

c  c 


where: 


N;  =  unknown  integer  ambiguity, 
f  s  measurement  frequency 

The  measured  beat  phase  based  pseudorange  in  meters  is: 

"  (<t>u,  +  Ni )  J  =  Pi  +  b«c  -  bGPSi  c  -  Apion 
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Pli,  =  -  (^li,  +  N, ) j  =  ^  +  buc  -  bGPSj  c  -  Apion  (4.8) 

J=Ni"Pi-buC  +  bGPS,C  +  APion  (4-9) 

Now  a  linear  combination  of  equations  (4.6)  and  (4.7)  is  made  to  remove  the 
change  in  range  due  to  the  ionosphere. 


Ppi,  +PLlj 
2 


Ppi,  (^Llj 
2 


(4.10) 


Ppi*  (^li,  +  Ni )  f  p.  +  buc  -  bGPSc  +  Apion  +  +  bnc  -  bGPSc  -  Api0 


(4.11) 


=PP‘‘  f  =P.+b,c-W  +  N|^-  (4.12) 

Now  the  ionospheric  free  DRVID  measurement  is  given  by 

Ppi.  -+h  f 

Y„= - j— £•  (4.13) 

The  ionospheric  free  computed  measurement  is  given  by: 

Computed  =  Pi  +  buC  —  bGPSi C  +  —  (4. 14) 

The  difficulty  now  becomes  solving  for  the  integer  ambiguity  Nj.  Nj  can 
be  estimated  as  part  of  the  state,  however,  care  needs  to  be  taken  to  account  for 
cycle  slips  on  the  phase  measurement. 
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4.2.1  Estimation  of  the  Integer  Ambiguity  Nj 

If  the  integer  ambiguity  is  to  be  estimated,  the  new  GEODE  state  would 
become: 

X  =  [x  y  z  x  y  z  CD  CR  b  b  N,  •••  Nj  •••  (4.15) 

where  n  is  the  number  of  GPS  satellites  providing  pseudorange  and  phase 
measurements  to  the  user  spacecraft  at  a  given  epoch.  The  observation  state 
relationship  is  defined  by: 

Y1l=G,[X(tl),tk]+el  (4.16) 

Since  the  observation  state  relationship  is  non-linear  we  take  the  partial  derivative 
of  the  computed  observation  equation  with  respect  to  the  state,  i.e., 

(with  b  =  buc  ) 

G,[X(t1),tl]  =  p,„^  =P+b-bGpsc  +  N,JL  (4.17) 

where 

P  =  \iXu~XGPS i)  +(yu-yGPSl)  (Zu  ~ ZGPSi  )  (4.18) 

Iik  _aGi[x(tk),tfc] 
j  ax(tk) 


(4.19) 
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The  difficulties  involved  in  this  implementation  are:  dynamically 
adjusting  the  size  of  the  state  to  reflect  the  number  of  GPS  spacecraft  being 
tracked  and  keeping  track  of  cycle  slips  to  reinitialize  the  appropriate  member  of 
the  state  when  a  cycle  slip  occurs.  Another  problem  with  this  implementation  is 
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the  added  computational  burden  and  complexity  of  adding  up  to  12  constants  to 
the  state,  depending  on  the  number  of  channels  of  the  GPS  receiver. 

4.2.2  Calculation  of  Apion  Using  Appl  and  ApL1 

Another  method  of  calculating  the  change  in  range  due  to  ionospheric 
effects  is  to  use  beat  phase  and  pseudorange  at  different  epochs.  In  this  manner 
the  change  in  the  ionospheric  effects  can  be  calculated  from  one  epoch  to  the 
next.  Differencing  equation  (4.9)  at  two  epochs  yields  [31]: 

Api,,  =(K, <4-21> 

where  0  is  the  initial  tracking  epoch  and  k  is  the  epoch  of  interest. 

Or 

APli.  =  -(Njt  -N° )-(p-  -p?)-(t>uC-b“c)+  (bgps. c - bGPSi c ) + ( Api0Ili  -Apion, ) 
Differencing  equation  (4.6)  at  the  two  epochs  yields: 

Ap‘t  =  (pf  -p')+(b:c-bIc)-(b‘psc-bJPSic)+(apL,i  -Apjj  (4.22) 
Assuming  the  integer  ambiguity  is  constant  between  epochs  and  adding  these  two 
“delta  range”  measurements  together  yields: 

ApPI,  -  Appl.  =  2  (Apfoni  -  Apfon.  )  (4.23) 

There  are  several  ways  to  account  for  Ap°n. .  One  is  to  model  the 
ionosphere  so  whenever  the  user  receiver  starts  tracking  a  new  GPS  spacecraft  the 
ionospheric ’’bias”  (Ap°011i )  is  calculated.  Then  with  Ap°n.  known: 
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ApL,  = 


App,,  ~  APl,. 


+  Apf0 


ion* 


(4.24) 


or 


Apfoni  =' 


(Pr^-PrO-fc-KO 


-+APL, 


(4.25) 


4.2.2.1  Model  Bias  DRVID 

In  Kubitschek  [61]  estimation  of  the  Ap°n.  bias  through  ionospheric 

modeling  using  GIM  maps  applied  in  conjunction  with  DRVID  provided 
excellent  ionospheric  error  estimation  compared  to  the  use  of  dual  frequency 
measurements.  Results  from  [61]  for  2  Feb  1997  GPS/MET  data  are  presented  in 
the  first  column  of  Table  4.1.  Estimation  of  the  Ap°n.  bias  using  the  GIM  maps 

requires  a  383  KB  data  file  that  is  updated  daily.  Uplinking  a  383  KB  file  to  a 
spacecraft  every  day  is  probably  unreasonable  overhead.  Another  approach  that 
proved  just  as  accurate  as  modeling  the  Apfon.  bias  in  [61]  is  not  estimating  the 
bias  but  allowing  the  OD  filter  to  estimate  it  along  with  the  GPS  satellite  clock 
bias.  Allowing  the  filter  to  estimate  the  Ap°n.  bias  works  because  the  bias  is 

constant  for  each  GPS  satellite  during  a  continuous  tracking  arc.  Unfortunately, 
real-time  filters  do  not  estimate  the  GPS  satellite  clock  biases  so  a  different 


approach  is  needed. 


Table  4.1  -  Error  Between  DRVID  Methods  and  Dual  Frequency 


3  Feb  1997 

RMS  Error  between 

dual  frequency  and 

bias  estimation  from 

bias  =  0 

PRN# 

modeling  (m) 

(m) 

PRN  1 

1.668 

1.876 

PRN  2 

1.888 

2.614 

PRN  3 

3.213 

3.297 

PRN  4 

1.635 

2.389 

PRN  5 

2.187 

2.651 

PRN  6 

1.121 

1.903 

PRN  7 

1.861 

2.035 

PRN  9 

1.22 

1.916 

PRN10 

1.985 

2.619 

PRN14 

2.305 

4.047 

PRN  15 

1.917 

2.463 

PRN  16 

2.576 

5.420 

PRN  17 

2.059 

2.264 

PRN  18 

2.515 

3.937 

PRN19 

1.672 

2.154 

PRN21 

1.927 

2.450 

PRN22 

5.74 

7.290 

PRN23 

1.752 

1.820 

PRN24 

2.138 

2.100 

PRN25 

3.171 

3.369 

PRN26 

0.996 

1.604 

PRN27 

1.235 

5.769 

PRN29 

1.009 

1.765 

PRN30 

0.944 

1.391 

PRN31 

1.333 

1.725 

Mean 

2.002 

2.835 

4.2.2.2  Zero  Bias  DRVID 

Another  less  computationally  burdensome  method  presented  itself  when 
plotting  the  measurement  elevation  angles  of  GPS/MET  with  respect  to  GPS 
satellites.  Since  GPS/MET  is  an  occultation  experiment  its  GPS  antenna  is 
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pointing  in  the  anti-velocity  direction  providing  reception  of  “setting” 
occultations.  See  Figure  2.2. 

Examination  of  the  elevation  of  the  GPS/MET  GPS  observations  indicates 
GPS  spacecraft  are  nominally  acquired  at  or  near  their  highest  elevation.  Figure 
4. 1  through  Figure  4.5  show  the  elevation  of  the  measurements  from  GPS 
satellites  PRN01,  PRN16,  PRN14,  PRN29  and  PRN10  taken  by  GPS/MET  on  4 
Feb  1997.  The  lighter  points  on  the  plots  distinguish  periods  when  the  OrbView- 
1  satellite  is  in  sunlight  (Day  in  the  legend  of  the  plots).  The  darker  points  are 
when  the  satellite  is  in  the  darkness  (Night). 

The  elevations  for  PRNs  1  and  29  show  they  are  nominally  acquired  at  or 
near  their  highest  elevation.  The  elevations  for  PRNs  16, 10  and  14  show  an 
elevation  increase  and  decrease  every  third  time  they  are  acquired.  Seven  out  of 
the  twenty-seven  GPS  satellites  tracked  by  GPS/MET  on  4  Feb  1997  exhibit  the 
elevation  decrease,  increase,  decrease  phenomenon  at  approximately 
-15°  elevation.  This  phenomenon  is  most  likely  due  to  occultation  of  the  GPS 


signal. 


Time  (hours) 

Figure  4.5  -  PRN10  GPS/MET  Elevation  Plot 


Most  of  the  GPS  tracking  arcs  start  with  GPS  satellites  at  their  highest 
elevation  with  respect  to  GPS/MET.  Therefore,  the  assumption  that  Ap°n,  ~0  is 

possible  since  the  change  in  range  due  to  ionospheric  effects  is  lowest  when  the 
GPS  spacecraft  are  at  their  highest  elevation.  Figure  4.6  shows  a  plot  of  the  dual 
frequency  ionospheric  correction  at  the  start  of  each  tracking  arc.  With  zero  bias 
DRVID  the  “bias”  at  the  start  of  each  tracking  arc  is  assumed  to  be  zero. 
However,  Figure  4.6  shows  the  mean  of  the  bias  at  the  start  of  each  tracking  arc  is 
2.12  m  with  several  significant  outliers.  The  two  largest  outliers  are  with  PRNs 
29  and  16.  Figure  4.3  shows  the  elevation  plot  for  PRN29  and  Figure  4.2  shows 
the  elevation  plot  for  PRN16.  As  seen  in  the  figures,  the  outliers  are  caused  by 
breaks  in  the  tracking  arcs.  At  the  time  of  each  of  these  tracking  arc  breaks  a 
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cycle  slip  occurs  on  the  phase  measurement  introducing  a  new  integer  ambiguity. 
Cycle  slips  also  disrupt  the  DRVID  calculations  since  the  phase  measurement 
does  not  have  the  same  integer  ambiguity  from  one  epoch  to  the  next.  The 
outliers  for  PRNs  10  and  14  are  similarly  explained.  Also  see  Figure  4.4  and 
Figure  4.5. 


Tracking  Arc  Start  Correction  Mean  =  2.12  m 
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Figure  4.6  -  Dual  Frequency  Ionospheric  Correction  at  GPS  Tracking  Arc  Start 


If  Apfon.  «  0  is  assumed,  equation  (4.25)  becomes. 


ApL,  = - « - - 


(4.26) 


Now  the  change  in  range  due  to  ionospheric  effects  can  be  approximated  without 
the  additional  computational  burden  of  estimating  integer  ambiguities  in  the  state 


or  by  estimating  the  initial  Ap°n,  bias.  To  apply  this  model  to  spacecraft  other 

than  GPS/MET,  application  of  the  DRVED  ionospheric  correction  should  only 
take  place  when  the  elevation  from  the  user  satellite  to  the  GPS  satellite  begins 
decreasing  and  when  the  elevation  decrease  begins  at  high  elevation  angles. 
Results  of  estimating  Apfon,  using  DRVED  on  the  GPS/MET  data  are 

shown  in  the  last  column  of  Table  4.1.  The  zero  bias  DRVID  estimation 
technique  suffered  only  a  0.833  meter  error  compared  to  the  model  bias 
technique.  Figure  4.7  shows  a  plot  of  the  largest  zero  bias  DRVID  errors 
compared  to  dual  frequency  shown  in  Table  4.1  and  Figure  4.8  shows  the  smallest 
zero  bias  DRVID  errors  shown  in  Table  4.1. 


!* 

■g  20 
5  0 

03  -20  H 
-40 


fa  ♦  GPS/MET  PRN22  1  ~  1  4 

i  h  \  \  •'  \  \  V  \  U  \ 


Figure  4.7  -  PRN  22  Zero  Bias  DRVED,  Dual  Frequency  Comparison 
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Figure  4.8  -  PRN  23  Zero  Bias  DRVTD,  Dual  Frequency  Comparison 


4.3  DRVID  in  GEODE 

The  implementation  of  DRVID  in  GEODE  requires  significant 
bookkeeping  of  the  GPS/MET  to  GPS  satellite  elevation  angles.  To  ensure  the 
measurement  epoch  with  the  smallest  ionospheric  error  is  chosen  to  provide 
p°j  and  (J)®!  for  equation  (4.26),  the  elevation  angle  and  the  phase  and 

pseudorange  measurements  from  the  previous  measurement  epoch  must  be 
available  at  the  current  epoch.  Then,  the  previous  elevation  can  be  compared  to 
the  current  elevation  and  when  the  current  elevation  is  smaller  than  the  previous 
elevation,  Ppj  and  (f^  are  recorded.  Then,  the  previous  epoch  is  the  zero  bias 
point.  Care  must  also  me  taken  to  ensure  cycle  slips  in  the  phase  measurements 


are  detected  during  any  given  tracking  arc.  If  a  cycle  slip  occurs,  the  zero  bias 
point  must  be  reset. 

Figure  4.9  shows  an  example  tracking  case  where  the  GPS  satellite 
elevation  with  respect  to  the  user  satellite  starts  at  85°  at  time  1,  then  increases  to 
90°  at  time  2  and  back  to  85°  at  time  3.  In  this  case,  the  pseudorange  and  phase 
measurements  from  time  2  will  be  used  as  p£j.  and  4)°,. .  The  measurements  at 

time  1  and  time  2  will  be  treated  as  a  single  frequency.  DRVID  will  be  applied  to 
the  measurement  at  time  3. 


Figure  4.10  shows  a  plot  of  DRVID  ionospheric  correction  estimates 
compared  to  dual  frequency  corrections  in  GEODE.  DRVID  appropriately 


estimated  the  necessary  ionospheric  corrections  in  most  cases.  PRNs  5, 27  and 
16,  cause  the  outliers  near  hours  11, 14  and  17  respectively.  Figure  4.12  through 
Figure  4.14  show  the  relationship  between  the  elevation  angle  and  the  error 
between  the  dual  frequency  ionospheric  correction  and  the  DRVID  ionospheric 
correction.  In  each  case  a  break  in  the  tracking  arc  and  cycle  slip  occurs  so  the 
zero  point  is  reset.  Since  the  breaks  occur  at  low  elevations,  single  frequency 
measurements  are  used  for  the  remainder  of  the  tracking  arcs.  The  vertical  lines 
on  Figure  4.12  through  Figure  4.14  indicate  when  the  satellite  is  in  sunlight  and 
darkness.  The  narrow  regions  are  darkness. 
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Figure  4.10-4  Feb  1997  GPS/MET  DRVID  Dual  Frequency  Comparison 
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Figure  4.11  -  PRN01  GPS/MET  DRVID  in  GEODE 
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Figure  4.12  -  PRN16  GPS/MET  DRVID  in  GEODE 
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4.3.1  DRVID  in  GEODE  Results 


The  application  of  the  zero  bias  DRVID  technique  to  the  GPS/MET  data 
in  GEODE  leads  to  an  accuracy  improvement  close  to  the  improvement  realized 
using  dual  frequency  data.  Table  4.2  shows  a  comparison  between  the  single 


frequency,  dual  frequency  and  DRVID  results  for  the  4  Feb  1997  data.  DRVID 


improves  the  3D  RSS  error  by  4.2%  for  the  case  without  high  rate  clocks  and 
4.5%  with  high  rate  clocks.  Additional  GPS/MET  DRVID  results  are  presented 
in  Chapters  5  and  6. 


Table  4.2  -  GEODE  With  Ionospheric  Correction  for  GPS/MET 


GPS/MET  -  4  Feb  1997 

Error  Mean  (m) 

RMS  Error  (ra) 

R 

I 

c  vy: 

R 

I 

c 

3D 

Without  High  Rate  GPS  Clocks 

Single  Frequency 

0.55 

-0.99 

-0.10 

2.74 

10.24 

2.15 

10.81 

Dual  Frequency 

0.11 

2.04 

-0.11 

2.48 

9.68 

2.12 

10.21 

DRVID 

0.24 

-1.14 

-0.11 

2.63 

9.72 

2.23 

10.32 

With  High  Rate  GPS  Clocks 

Single  Frequency 

0.28 

-3.59 

-0.04 

2.22 

7.23 

0.77 

7.60 

Dual  Frequency 

-0.08 

-1.38 

-0.04 

2.10 

6.14 

0.86 

6.55 

DRVID 

-0.03 

-3.54 

-0.04 

1.97 

6.95 

0.90 

7.28 

Zero  bias  DRVID  ionospheric  corrections  are  also  applied  to  the  SA  free 
T/P  data.  Table  4.3  shows  the  3  cm  3D  RSS  position  error  improvement.  The 
reason  the  accuracy  improvement  is  proportionally  smaller  (only  2.5%)  than  the 
improvement  realized  on  the  GPS/MET  data  is  the  ionospheric  effects  on  the  T/P 
data  are  smaller  due  to  T/P’s  higher  orbit. 


Table  4.3  -  GEODE  With  Ionospheric  Correction  for  T/P 


T/P -5  May  2000 

Error  Mean  (m) 

RMS  Error  (m) 

R 

I 

C 

R 

I 

C 

3D 

Single  Frequency 

0.21 

-0.04 

0.01 

0.35 

1.06 

0.44 

1.20 

DRVID 

0.20 

0.01 

0.01 

0.35 

1.02 

0.45 

1.17 

While  the  error  RMS  and  RSS  improvements  are  near  those  of  using  dual 
frequency  ionospheric  corrections,  a  comparison  of  the  execution  time  required 
by  GEODE  with  and  without  DRVID  shows  a  12.1%  computational  burden 
increase.  This  significant  increase  is  a  result  of  the  need  to  calculate  and  store  the 
elevation  angle  for  each  satellite  being  tracked  at  every  measurement  epoch. 
Calculating  the  GPS  elevation  angle  requires  computation  of  the  position  of  each 
GPS  satellite  in  view,  therefore,  the  GPS  satellite  ephemerides  must  be  evaluated 
for  each  elevation  calculation.  The  computational  burden  is  reduced  significantly 
if  all  satellites  in  view  are  processed  for  the  measurement  update  (ALL)  instead  of 
cyclically. 


4.4  Summary 

Differenced  Range  Versus  Integrated  Doppler  (DRVID)  [29]  is  a 
technique  utilizing  the  difference  in  the  way  the  pseudorange  and  phase 
measurements  are  affected  by  the  ionosphere.  Since  the  ionosphere  advances 
phase  and  delays  range  measurements,  a  linear  combination  of  both 
measurements  can  remove  ionospheric  errors  to  first  order. 

There  are  three  methods  by  which  DRVID  can  be  applied  to  estimate  GPS 
pseudorange  ionospheric  corrections.  First,  the  phase  integer  ambiguity  can  be 


estimated  as  part  of  the  state.  Integer  ambiguity  resolution  is  not  suited  to  real¬ 
time  systems  and  therefore  is  not  considered  here.  The  second  application  of 
DRVID  involves  estimation  of  an  initial  ionospheric  bias  correction  using  a 
model  such  as  JPL’s  GIM.  Model  bias  estimation  shows  promise,  but  the 
overhead  involved  in  generating  the  bias  estimates  and  updating  the  model  is 
probably  unreasonable  for  a  real-time  system.  Application  of  the  zero  bias 
DRVID  technique  to  estimate  the  GPS  pseudorange  ionospheric  correction  results 
in  an  accuracy/precision  improvement  over  using  single  frequency  measurements 
with  no  ionospheric  correction.  The  zero  bias  DRVID  technique  provides 
accuracy  very  near  dual  frequency  correction  but  with  a  12.1%  increase  in 
computational  burden. 


CHAPTER  5 


Dynamic  Model  Compensation  (DMC) 

The  dynamic  models  used  to  propagate  satellite  ephemeris  are  always 
approximations  of  the  actual  forces  acting  on  the  satellite.  No  matter  how 
detailed  and  precise,  models  will  always  fall  short  of  describing  the  real 
dynamics.  In  real-time  satellite  OD  the  dynamic  model  is  used  to  propagate  the 
satellite’s  position  and  velocity  between  measurement  epochs.  The  dynamic 
model  might  also  be  used  to  propagate  the  satellite’s  ephemeris  for  prediction 
purposes.  If  the  dynamic  model  does  not  predict  the  satellite’s  motion 
accurately  enough,  the  filter  in  an  OD  scheme  can  diverge.  Care  needs  to  be 
taken,  in  a  real-time  OD  system  to  balance  the  accuracy  requirements  of  the  time 
update  of  the  state  with  the  overhead  of  the  propagator.  In  choosing  a  finite 
approximation  of  the  satellite  dynamics  there  will  always  be  forces  in  the 
dynamic  model  that  will  be  unmodeled  or  mismodeled. 

Dynamic  Model  Compensation  (DMC)  is  a  process  noise  formulation  that 
assumes  a  dynamical  system  is  subject  to  accelerations  not  included  in  the 
systems  dynamic  model  and  which  possess  a  random  element  [9].  These 
accelerations  have  been  called  “fictitious”  [62],  “augmenting”,  “compensative” 
[9]  and  “empirical”[4].  Regardless  of  what  they  are  called,  they  are  estimated  in 
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the  filter  and  their  formulation  includes  the  development  of  process  noise 
covariance  values.  DMC  was  used  as  early  as  the  1970’s  [63]  but  the 
formulation  shown  here  follows  the  derivation  in  Cruickshank  [64]  and  [9]. 

In  this  chapter,  DMC  is  first  introduced  through  a  simple  one-dimensional 
example.  Then  the  development  is  extended  to  the  three-dimensional  satellite 
OD  problem  where  accelerations  are  estimated  in  Cartesian  XYZ  coordinates. 

Next,  DMC  is  refined  and  accelerations  are  estimated  in  Radial,  In-track  and 
Cross-track  (RIC)  coordinates.  Finally,  XYZ  and  RIC  DMC  are  employed  in 
two  simulations  and  in  GEODE  and  accuracy/precision  and  computational 
burden  results  are  presented. 


5.1  Dynamic  Model  Compensation  (DMC) 

Suppose  the  particle  in  Figure  5.1  is  moving  with  constant  velocity  and  the 
particle’s  position  is  to  be  estimated  from  range  observations. 

Vo 

• - ► 

- 1 - 

Xo 

Figure  5.1  -  Particle  Moving  at  Constant  Velocity 


The  dynamic  model  for  this  system  can  be  expressed  in  state  space  form  as: 


(5.1) 


leading  to  the  derivative  in  state  space  form. 
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"0 

r 

X 

0 

0 

X 

X  =  AX  = 

The  state  transition  matrix  for  this  simple  system  is  formed  from 

6(t,t0)  =  A<D(t,t0) 


(5.2) 


(5.3) 


where  the  initial  condition  is  4>(t0,t0)  =  I.  A  method  of  calculating  the  state 
transition  matrix  for  this  simple  system  is  Laplace  transforms.  Here 

<t(t,t0)  =  L-,[(sI-A)-1]  (5.4) 

The  state  transition  matrix  is  therefore 


<t>(t,t0) 


1  (t-‘o) 

0  1 


(5.5) 


If  stochastic  forces  perturb  the  movement  of  the  particle,  either  random  (possibly 
time  correlated)  and/or  deterministic,  estimation  under  the  constant  velocity 
assumption  is  not  optimal.  DMC  can  be  used  to  improve  estimation  performance 
in  such  a  situation.  Assume  that  the  unmodeled  acceleration  co(t)  can  be  modeled 
as  a  first  order  Gauss-Markov  process 

cb(t)  =  pco(t)  +  u(t)  (5.6) 


where: 


co(t)  s  compensative  acceleration 

u(t)  =  uncorrelated,  stationary  Gaussian  process  (white  noise)  with 
E[u]  =  0  and  E[uTu]  =  a2u  (standard  deviation  of  the  forcing  noise) 

p  =  —  where  x  is  the  correlation  time  and  is  considered  constant 
x 


Now  the  state  can  be  augmented  with  the  compensative  acceleration 
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X  = 


x 

x 

(0 


(5.7) 


and 


X 

0 

l 

o' 

X 

'o' 

x  = 

CO 

= 

0 

0 

1 

X 

+ 

0 

-Pco+u(t) 

0 

0 

-p_ 

(0 

l 

u(t)=AX+Bu(t) 


(5.8) 


The  general  solution,  found  using  convolution,  has  both  deterministic  and  random 
components 


X(t)  =  <Dw(t,t0)X(t0)  +  J‘  «Dw(t,t0)Bu(t)dt 


(5.9) 


It  1 \ 

Deterministic  Random 

The  state  transition  matrix  Ow(t,t0)  can  again  be  found  from  <j>(t,t0)  =  Ad>(t,t0). 
The  Laplace  transform  method  can  again  be  used,  but  inverting  a  3x3  or  larger 
matrix  can  be  a  challenge.  The  Symbolic  Toolbox  in  Matlab  or  Mathematica  can 
easily  be  used  for  larger  matrices. 


<MU0)  = 


1  At  TAt  +  x2(e-At/T-l) 
0  1  x(l-e~At/T) 

0  0 


-At/T 


(5.10) 


Assuming  the  initial  stochastic  state  is 


X(t0)  = 


0 

0 


L°V 


(5.11) 


the  deterministic  component  of  equation  (5.9)  is 


Xdel  —  <I)w(t,t0)X(t0)  — 


co0xAt  +  co0x2  (e  At/T  - 1) 
©0x(l-e-At/T) 


o)0e 


-At/t 


116 

(5.12) 


Now  the  state  propagation  model  is  also  augmented  and  the  deterministic  portion 
of  (5.9)  is  included. 


'x(t)' 

to0xAt  +  co0x2(e_Al/T  - 1)  +  v0At  +  x0 

X(t)  = 

x(t) 

r= 

co0x(l-e-A,/T)  +  v0 

_ro(t)_ 

co0e-A,/T 

The  process  noise  covariance  can  be  formulated  from  the  random  terms  in  (5.9), 
i.e.  f '  <tw  (t,t0)Bu(t)dt .  The  calculation  of  the  process  noise  covariance  matrix 

Jto 

over  the  time  interval  At  =  t  - 10  is: 

Qw  =  J‘  <Dw(t,x)BE[u(t)uT(x)]BT0)^t,x)dx  (5.14) 

or  since  E[u(t)uT(x)]  =  a2 

Qw  =  f  <Dw(t,x)Ba2BTC(t,x)dx  (5.15) 

J‘o 

where  x  is  the  time  integration  variable,  not  the  stochastic  process  correlation 
time.  The  resulting  integral  is: 


"x(t-T)+x2  (e~(tT)/T  —  lj 

x(t-T)  +  x2(e“(tT)/T-l) 

Q.  = 

x(e-(t'T)/T) 

x(e-(,'T)/t) 

% 

l0 

e-(t-T)/x 

e-(t-T)/x 

(5.16) 


where  T  is  substituted  for  the  time  integration  variable  x. 
Evaluation  of  the  integral  in  equation  (5.16)  results  in: 


where: 


Q 

Q 

Q 


W(l,l) 

Qw(1,2) 

Qw(1,3) 

w(2,l) 

Qw(2,2) 

Qw(2,3) 

w(3,l) 

Qw(3,2) 

Qw(3,3) 
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(5.17) 


Addition  of  the  process  noise  covariance  to  an  Extended  Kalman  Filter 
(EKF)  is  rudimentary.  The  algorithm  is  unchanged  except  for  the  time  update  of 
the  state  error  covariance  Pk  [Bom,  2000  #252].  The  new  formulation  for  Pk  is: 

I*k+1  =  ®(^k+l’^k  (tk+l’tk)“*“Qw  (5-19) 

At  =  tk+,  -  tk  in  the  formulation  of  Qw. 

The  application  of  DMC  to  a  real-time  satellite  OD  system  using  GPS  is 
more  complex  since  the  state  includes  constants  from  the  dynamic  model,  receiver 
clock  bias  and  receiver  clock  drift  terms.  Therefore,  the  state  transition  matrix 
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and  process  noise  covariance  matrices  explained  above  are  updated  in  section  5.2 
to  include  the  additional  terms. 

5.2  XYZDMC 

The  XYZ  DMC  derivation  shown  below  follows  closely  the  derivation  in 
Cruickshank  [9].  As  shown  in  equation  (5.6)  the  stochastic  accelerations  are 

modeled  dynamically  by  a  Langevin  equation  [9].  Now  3  in  equation  (5.6) 
becomes: 

~1/t  0  0" 

p=  0  1/t  0  (5.20) 

0  0  1/t 

The  Gaussian  process,  u(t),  is  uncorrelated  in  time  with  zero  mean  but  now  has 
constant  variance  of: 

0  0 ' 

qu=  0  0  (5.21) 

_°  0  oJ_ 

P  and  qu  are  assumed  diagonal  for  convenience  only,  indicating  the  accelerations 
in  the  three  directions  are  uncorrelated  with  each  other,  t  and  gu  are  also 
assumed  equal  in  each  of  the  three  directions  [9]. 

The  state  vector  is  augmented  with  the  compensative  accelerations,  cox , 
coy ,  ooz  as  shown  in  equation  (5.22).  The  deterministic  components  augment  the 
time  propagation  of  the  state  as  in  equation  (5.23). 


The  original  New  state  vector 

state  vector  in 
GEODE 


X 


y 

z 

X 


(5.22) 


co„ 


CO, 


z  Jl3xl 
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New  state  augmented  with  deterministic  portion  of 
compensative  accelerations 


X  +  G)XT  At  +  C0xT2  (e_At/T 
y  +  coyT  At  +  coyT2  (e_At/x 
z  +  oozx  At  +  cozx2  (e~At/x 
x  +  coxx  (l-e_At/T  ) 
y  +  CDyT  (l-e~At/x  ) 
z  +  cdzx  (l-e_At/x  ) 

CD 

CR 

b 

b 


-O' 

-0 

-0 


or  e 

xk-i 


-At/x 


<Ve 


-At/x 


coz  e 

zk-i 


-At/x 


(5.23) 

13x1 


5.2.1  XYZ  DMC  State  Transition  Matrix 

The  state  transition  matrix  is  also  augmented  to  include  the  compensative 
acceleration  terms. 

(Jj  _  [^xliOxlO 

L  [°Lo 

The  state  transition  matrix  of  the  original  state  is  shown  in  equation  (5.25)  and  is 
developed  in  Lee  [37]. 


[^coRV  ] 

[i'.l 


10x3 

3x3 


<D  = 


ax 

axn 


(5.24) 
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<t>  = 


’  5R(tk)  " 

’  5R(tk) " 

3R(tk) 

3R(tk) 

_3R(tk.,)_ 
’  SR(tk) 

3x3 

.«R(tk.,). 
’  5R(tk)  " 

3x3 

-^ACD(tk,)_ 

r  5R(tk>  i 

3x1 

-5ACR(tk.,)I 

2R(tk) 

.3R(tk„). 

3x3 

.3R(tk_,)_ 

3x3 

1 

o 

% 

_ 1 

3x1 

aAca  ,) 

L.  R  x  k-1  '  — *3x 

[oL 

[0]M 

[°L 


[o]„ 

[o]„ 

K 


[oL 

(5.25) 


dAC„(tt) 

3ACR(tk.,)J 

[oL 


[oL, 

[OL 

[oL 

[0]., 

3Ab(tt) 

.SAb(tw)J 


The  state  transition  matrix  of  the  compensative  accelerations  with  respect  to 
themselves: 


0>  = 


dcox 

d(ox 

do\ 

s 

X 

o 

5COVo 

d(°c0 

8coy 

dcoy 

dcoy 

dcox 

x0 

5®c0 

8(oz 

8(oz 

dcoz 

do\, 

do\ 

8(or 

M> 

-At/t 

0 

0 


0 

j-ai h 

0 


0 

0 

s-At/x 


At  —  tk+i  tk 


(5.26) 


Again, 


The  state  transition  matrix  of  the  compensative  accelerations  with  respect  to  the 
original  state  is: 
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(5.27) 


(5.28) 


and 


dx  6x  dx 

5®x0  dCOy0  eCOZo 

fo  i=  ^  _  3_  JL.  JL 

0)V  3co0  6<ox  6(0  6(0 

U  x0  y0  z0 

6i  6z  6z 

6ax  6(ov  d(o7 

L  xo  yo  zo 


(5.29) 


Therefore, 
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Atx  +  x2(e'At/T  -l)  0  0 

0  Atx  +  x2(e~At/T  -l)  0 

0  0  Atx  +x2(e-At/T  -l) 

x  (l-e~At/T )  0  0 

4>„rv=  0  t  (l-e-^'-  )  0 

0  O  x  (l-e-A,/T ) 

[OOO" 

0  0  0 
0  0  0 
[O  0  0 


5.2.2  XYZ  DMC  Process  Noise  Covariance 

The  process  noise  covariance  matrix,  Qxyz,  is  propagated  over  the  time 
interval  to  to  t  through  evaluation  of  the  integral  in  equation  (5.31). 

"o2  0  0 1 

Qxyz =  o  a2  0  f‘d>(t,T)BBTd)T(t,T)dT  (5.31) 

Jt° 

L°  °  °i\ 

Here,  B  is  a  matrix  of  the  form: 

"0  0  0  0  0  0  0  0  0  0  1  0  olT 

B=0000000000010  (5.32) 

0000000000001 


au2  0  0 

S=  0  a2  0 

o  0  at 


(5.33) 
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Then  the  DMC  contribution  to  Qxyz  becomes: 


]  [^©R  ]3x3 

[°L 

[k-IKL 

sf 

[®tov][®coR  ]3x3 

[®wv][®coV  ]3x3 

[°L. 

[®.v][®.£, 

Jt0 

[°L 

[°L, 

[°L. 

[°L, 

K1K.L 

Now  let 


Qxyz  (1,1) 

Qxyz(l,2) 

[°L 

Qxyz(l,4) 

n 

Qxyz(2,l) 

Qxyz(2,2) 

[°L 

Qxyz(2,4) 

VXyz  ~ 

[°L.3 

[®L 

Qother4x4 

[°L 

Qxyz(4,l) 

Qxyz  (4,2) 

[°L 

Qxyz(4,4) 

(5.34) 


where  the  (3,3)  sub-matrix  Qothen  contains  process  noise  contributions  for  Cd,  Cr, 
b  and  b .  Also, 


Qxyz(l,l)  [^oRH^oR^xjdT 


(5.35) 


The  rest  of  Qxyz  is  defined  as  in  equation  (5.35).  Now 

[d>wR][G>*R]  =  [Atx  +  x2(e'At/T  l)]  [ Atx  +  x2(e~At/T  -l)]l 

and  the  integral  in  equation  (5.35)  is  evaluated  yielding 


3x3 


-xyz(l.l) 


(5.36) 


— t5  (l-e'2At/T  )+x4At(l-e"A,/t  ) -  x3 At2  + — x2 At3 1 S  (5.37) 


The  other  sub-matrices  of  equation  (5.34)  are  determined  in  like  fashion. 

Q-HU)  )+W 


(5.38) 


Qxyz(l,2)  Qxyz(2,l) 


Qx^(..-)  =  H Ate'"'’ 


-xyz(2,2) 


■-+ 2e-A,/T  --e'2A,/T 

v  2  2  / 


+  x2At 


Qxyz(2,4)  -Qxyz(4,2)  -x2fT(1  +  e  ^  )~e  ^ 


The  development  of  the  Qother  matrix  can  be  found  in  [37], 


Qother 


crLgAt  0 

0  ^LAt 


0  0 

0  0 


0  0 

0  0 

a-At3  ,  CT2At2 

- +  CT2At  - 

o  D  ^ 


a2  At2 
0 


a2  At 

b 
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(5.39) 

(5.40) 

(5.41) 

(5.42) 


(5.43) 


5.3  Radial,  In-track,  Cross-track  (RIC)  DMC 

The  development  that  follows  is  motivated  by  the  desire  to  estimate  the 
compensative  accelerations  in  the  radial,  in-track,  and  cross-track  directions.  RIC 
DMC  also  provides  a  more  intuitive  approach  to  tuning  the  time  correlation 
coefficient  (x)  and  standard  deviation  of  the  forcing  noise  (au ) .  The  state  is, 

once  again,  augmented  with  the  compensative  accelerations  as  shown  in  equation 
(5.44).  The  deterministic  portion  of  the  compensative  accelerations  now  must  be 
rotated  from  the  RIC  frame  to  the  XYZ  frame  to  be  added  to  the  state.  Equation 


(5.46)  shows  the  rotation  of  the  position  components  and  equation  (5.47)  shows 
the  rotation  of  the  velocity  components. 
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The  original  state 
vector  in  GEODE 


New  state  vector 


New  state  augmented  with  deterministic 
portion  of  compensative  accelerations 


X  = 


X 

x  +  Ax 

y 

y  +  Ay 

X 

z 

z  +  Az 

y 

X 

x  + Ax 

z 

y 

y  +  Ay 

X 

z 

z  +  Az 

y 

X  = 

cD 

(5  .44)  X  = 

cD 

z 

Cr 

Cr 

cD 

b 

b 

Cr 

b 

b 

b 

b  j 

10x1 

(0R 

®i 

coR  e-At/TR 

Kk-1 

<j0i  e“At/T‘ 

*k-l 

C0c 

13x1 

coc  e"A1/tc 

(5.45) 


Jl3xl 


Ax 

®RTR  At +  ®RTR  (e  R 

A*xyz  = 

Ay 

=  T 

ARIC-»xyz 

cOjTjAt  +  tOjTj  (e_At/T'  -l) 

Az 

tocrcAt  +  o)cTc(e'At/Tc  -l) 

(5.46) 


aL  = 


(■*  -At/tn  \ 

‘Ax' 

®R^(l-e  R) 

Ay 

=  T 

ARIC-»xyz 

Az 

®CTC  V  e  j 

(5.47) 


TRic_>Xyzis  the  transformation  matrix  from  RIC  coordinates  to  xyz  coordinates. 


The  equations  describing  TRic->.xyz  can  be  found  in  Lee  [37]. 
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5.3.1  RIC  DMC  State  Transition  Matrix 

The  [<DX]  and  [d>m]  matrices  are  the  same  for  the  XYZ  and  RIC  versions 

of  the  state  transition  matrix.  However,  the  state  transition  matrix  of  the 
compensative  accelerations  with  respect  to  the  original  state  must  be  rotated  to 
RIC  coordinates. 


Here, 


^(oRV  — 


dx_ 

da 

dx_ 

da 

gcD 

da 

dC* 

da 

db' 


3x3 


J3x3 


1x3 


1x3 


da 

db 

da 


-Ilx3 


Jlx3 


d t 

da 

[<>1 


[^coR  ]3x3 

]3x3 

.  [0Lx3 


(5.48) 


dr 

da 

dr__ 

da 

—  = 

da 


dr  d7mc 
c^jc  da 

(5.49) 

y  ^IC 

■*-RIC-»xyz 

(5.50) 

yT  ^,c 

RIC<-xyz  d- 

(5.51) 

with  a  similar  derivation 


127 


=tt 

—  ARIC<-xyz  ^ — 
0(0  0(0 

therefore  <f>C)RV  becomes, 


AtxR+4(e-^-l) 


TTRIC«-*yz 


0 

0 


0 

Atx,  +  T*(e"At/T'  -l) 
0 


0 

0 

Atxc+x^(e-^-l) 


®d>RV  — 


RIC<-xyz 


0 

0 


0 

x,(l-e-^) 

0 


0  0  0 
0  0  0 
0  0  0 
0  0  0 


and  in  simplified  form: 


Tt  d> 

ARIC<-xyz  ^  coRjyc 

1 

P4 

3 

e 

®coRV  ” 

tt 

A  RIC<-xyz  ^  coVrjc 

= 

^toV 

.  t°L  . 

J°L. 

5.3.2  RIC  DMC  Process  Noise  Covariance 

In  RIC  DMC  separate  standard  deviations  of  the  forcing  noise  (as)  are 

used  instead  of  equal  values  on  the  diagonal  as  in  XYZ  DMC.  Therefore,  the  S 
matrix  takes  the  form 


S  = 


0  0 
of  0 
0  g2c 


(5.52) 


(5.53) 


(5.54) 


(5.55) 
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causing  Qric  to  keep  the  form  of  equation  (5.56)  since  the  S  matrix  cannot  be 
factored  out. 


Qkic  -  L 


KM®.,} 

KJSK£ 
[°L 
KM®™] 


|T 

3x3 


13x3 


T 

3x3 


K]sKl 
KM®..] 
[°L 
KM®..] 


T 

3x3 

T 

3x3 


T 

3x3 


[0} 

[0] 

[0] 

[0] 


3x4 


3x4 


4x4 


3x4 


K]SKl 
K]SKJ 
[®l 
KM®.} 


T 

3x3 

T 

3x3 


T 

3x3 


UT  (5.56) 


With  the  new  formulation  of  the  state  transition  matrix,  the  upper  left  3x3  sub¬ 
matrix  of  Qric  becomes: 

[®or  ]3x3  =  [TRic^xyz^MR^  ]s[TRIC<_xy2<I)(()RRic  ] 

jSf^oR  ]3x3  =  [TRIC^xyz(t)(oRRic  ]s[(l)(oRR]cTRIC<_xyz  J 

[^ojrJS^^r]^  =[TRiC<-xyz^>fflRwcS<l>0)RRlcTRIC<_xyzJ  (5.57) 

Since  4>aR  (defined  in  equations  (5.53)  and  (5.54))  and  S  are  both  diagonal 

[^mR  ]3x3  =  [TRIC<_xyzS<broRR[c  (b0)RRICTR[C^xyz  J 

The  same  development  applies  to  the  other  members  of  the  Qric  matrix. 

With  Tric  =  TRIC<_xyz ,  Qricci.i)  becomes 

Qric(U)  =  J,0  [TRics4)(0RRIC4>fflRRlcTRIC]3x3dT  (5.58) 


Since  the  transformation  matrix  Tric  and  S  are  both  independent  of  the  variable  of 
integration  T, 


Q 


RIC(U) 


(5.59) 


Now  let 
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and  define  Q  as 


Q(u) 

Q(l,2) 

[°L 

Q(1.4) 

n  - 

Q(2,l) 

Q(2,2) 

[°k. 

Q(2,4) 

[°]4X3 

[<>]« 

[°L. 

[°]4X3 

_  Q(4,l) 

Q(4,2) 

[«L 

Q(4,4) 

Now  Qric  becomes, 


TRicQ(1,1)Tric 

TrXcQ(1,2)Tric 

[°L 

TricQ(1,4) 

n 

r^RIcQ(2,l)r^RIC 

TRicQ(2,2)r^RIC 

[°k< 

TricQ(2,4) 

Vric  — 

[°L, 

[°Lj 

Q  other 

[°L, 

Q(4,l)^RIC 

Q(4,2)Tric 

[ok. 

Q(4,4) 

(5.60) 


(5.61) 


The  upper  left  term  is  shown  as  an  example.  The  other  terms  are  formulated  in 
like  fashion. 


~[^tTx +Tx(e  A* ,x  l)][Atxx  +xx  (e  x  l) 


1  0  0 
0  1  0 
0  0  1 


(5.62) 


In  the  (1,1)  term  of  <b(tfRRIC(l)^RR[c  set  X  s  R,  in  the  (2,2)  term  X  =  I  and  in  the  (3,3) 
term  X  =  C.  Now, 

Q(i,i)x  =  £  [ Atxx  +  4  (e"At/Tx  - 1)]  [ AtTx  +  4  (e"At/Tx  -  l)]dT  = 

J'[4(e"2(t‘T)/Tx-2e-(,-T)/tx+l)+ 

2(t-T)4  (e"(,'T)/Tx  -l)+(t2  -2tT+T2)x2  ]dT  = 

(j4  (l-e"2At/Tx  )+  4At(l-e-At/Tx  )-x3xAt2  +  ^x2a/ 


(5.63) 


And  again, 
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Q  other 


^dragAt 

0 

0 

0 


0 

°LAt 


o 

o 


0 

0 


crAt  2 
“V-+abAt 


CrAt 

b 


0 

0 

a- At2 

b 


o?At 

b 


Finally, 


Tr,cQ(U)Tric 

TricQ(1,2)Tric 

[°L 

TricQ(1,4) 

TricQ(2,1)TrIC 

TricQ(2,2)Tric 

[°L 

TricQ(2,4) 

[ok. 

[°L 

Q  other 

[°L. 

Q(4,1)Tric 

Q(4,  2)Tric 

[°L 

Q(4,4) 

(5.70) 


(5.71) 


5.4  Single  Satellite  Range  and  Range  Rate  Simulation  (Simulation  1) 

Two  simulations  are  developed  in  Matlab  to  compare  the  results  of  XYZ 
DMC  and  RIC  DMC.  Both  simulations  consist  of  a  user  satellite  in  LEO 
(modeled  after  the  QuikScat  orbit).  Satellite  Toolkit  (STK)  is  used  to  generate  a 
“truth”  trajectory  using  a  70x70  gravity  model,  atmospheric  drag,  solar  pressure 
and  third-body  effects  from  the  Sun  and  Moon.  The  trajectory  is  propagated  over 
20,000  seconds.  Next,  tracking  satellites  in  geostationary  orbits,  with  equally 
spaced  longitudes,  are  used  to  generate  simulated  observations.  The  first 
simulation  uses  three  “TrakSats”,  each  separated  by  120°  of  longitude.  Range 
and  range  rate  measurements  are  taken  at  20-second  intervals  with  only  one 
TrakSat  measurement  at  each  observation  epoch.  Then,  zero  mean,  standard 
deviation  1  m  Gaussian  noise  is  added  to  the  range  measurements  and  zero  mean, 


standard  deviation  0.1  m/s  Gaussian  noise  is  added  to  the  range  rate 
measurements. 

An  Extended  Kalman  Filter  (EKF)  is  then  implemented  to  process  the 
range  and  range  rate  observations.  The  dynamic  model  used  by  the  EKF  consists 
of  two-body  and  J2  gravity  terms  and  an  exponential  drag  model.  The  state 
transition  matrix  is  integrated  along  with  the  position  and  velocity  and  includes 
two-body,  J2  and  exponential  drag  terms. 

5.4.1  Simulation  1  No  Process  Noise 

The  initial  EKF  implementation  did  not  include  DMC  nor  did  it  include 
any  form  of  process  noise.  Figure  5.2  shows  the  filter’s  estimated  trajectory 
compared  to  the  STK  truth  data.  The  standard  deviation  in  the  error  plot  is  the 
square  root  of  the  estimated  state  error  position  variance  term.  The  covariance 
approaches  zero  since  no  process  noise  is  used. 


0  1  2  3  4  5  6 


0  1  2  3  4  5  6 


Figure  5.2  -  Simulation  1  No  Process  Noise  Results 


5.4.2  Simulation  1  GEODE  Process  Noise 

Next  the  process  noise  formulation  used  in  GEODE  is  added  to  the 
simulation  filter,  t  and  the  a  s  are  tuned  to  minimize  the  3D  RSS  error  of  the 
filtered  solution  when  compared  to  the  truth  data.  The  error  and  standard 
deviation  statistics  are  shown  in  Figure  5.3.  Figure  5.3  shows  that  by  adding 
process  noise,  precision  and  accuracy  improved  and  the  variances  appear  more 
reasonable. 


3D  RSS  Error  =  75.39  m 


Figure  5.3  -  Simulation  1  GEODE  Process  Noise  Results 
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5.4.3  Simulation  1  XYZ  DMC 

Next  the  EKF  is  implemented  with  XYZ  DMC  and  t,  ctu,  c^,  on  and  ctcd 

are,  again,  tuned  to  produce  the  lowest  3D  RSS  error  compared  to  the  STK  truth 
trajectory.  Adding  the  estimation  of  the  compensative  accelerations  and  the  DMC 
process  noise  significantly  improved  the  filter’s  precision  and  accuracy.  Figure 
5.4  shows  the  error  and  standard  deviation  statistics  for  the  XYZ  DMC  case. 


3D  RSS  Error  =  68. 1 2  m 


Figure  5.4  -  Track  Simulation  1  XYZ  DMC  Results 
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5.4.4  Simulation  1  RIC  DMC 

The  final  step  in  simulation  1  is  to  convert  the  compensative  accelerations, 
state  transition  matrix  and  process  noise  covariance  matrix  to  RIC  coordinates  and 

tune  the  ts  and  as  (rR,  n,  Tc,  aR,  oh  oq,  ap,  an  and  ocd)-  Results  are  shown  in 
Figure  5.5. 


3D  RSS  Error  =  43.35  m 


L— 

o 


Figure  5.5  -  Simulation  1  RIC  DMC  Results 


Table  5.1  shows  a  summary  of  the  results  of  simulation  1.  Adding  process 
noise  to  the  filter’s  covariance  matrix  significantly  improves  its  accuracy  and 
precision.  Additional  accuracy  and  precision  are  gained  by  estimating 


compensative  accelerations  and  by  applying  RIC  DMC.  RIC  DMC  improved 
the  filter’s  error  RSS  by  36.4%  compared  to  XYZ  DMC. 


Table  5.1  -  Summary  of  Results  with  3  TrackSats  and  1  m  Noise  on  Range 


Simulation  1 

Error  Mean  (m) 

Error  RMS  (m) 

R 

I 

.  c 

R 

I 

C  ,  . 

3D 

No  Process  Noise 

-23.5 

92.1 

-22.3 

57.8 

141.6 

108.4 

186.0 

GEODE  Process  Noise 

-4.6 

7.5 

-20.5 

50.0 

40.5 

39.3 

75.4 

XYZ  DMC 

-0.7 

6.2 

-21.7 

38.2 

39.0 

40.7 

68.1 

RIC  DMC 

1.6 

2.7 

-12.5 

19.5 

19.9 

33.2 

43.3 

5.5  Multi-satellite  Range  Simulation  (Simulation  2) 

The  second  simulation  uses  the  same  STK  trajectory  as  truth.  However,  in 
this  simulation  four  TrakSats  are  used,  each  spaced  90°  apart  in  longitude.  Range 
measurements  are  then  calculated  artificially  from  all  four  satellites  at  20-second 
intervals.  In  simulation  2  there  are  four  range  observations  at  each  observation 
epoch  and  no  range  rate  observations.  Obviously,  measuring  range  from  all  four 
satellites  at  once  is  not  possible  due  to  interference  by  the  Earth.  The  additional 
satellite  is  added  and  measurements  from  all  satellites  at  once  are  used  to  improve 
the  simulation’s  observation  geometry.  The  improved  geometry  allows  a 
significant  improvement  in  accuracy  and  precision.  The  3D  RSS  error  without 
DMC  is  158  m.  When  RIC  DMC  is  applied  the  3D  RSS  error  is  3.1  m.  As  shown 
in  Table  5.1,  the  best  result  using  RIC  DMC  in  simulation  1  is  43.3  m,  3D.  Since 
the  objective  of  this  simulation  is  to  determine  if  transforming  DMC  to  the  RIC 
orientation  improves  solution  accuracy,  mean  of  zero,  standard  deviation  50  m 
Gaussian  noise  is  added  to  the  range  measurements  instead  of  1  m  noise.  The 


filter  result  without  process  noise  is  a  3D  RSS  error  of  169  m  but  with  GEODE’s 
process  noise  formulation  the  best  result  is  50.3  m  3D.  Figure  5.6  through  Figure 
5.9  show  error  plots  for  the  various  methods  and  Table  5.2  shows  a  summary  of 
the  error  statistics  for  the  various  methods. 


3D  RSS  Error  =  169.30  m 


Figure  5.6  -  Simulation  2  No  Process  Noise  Results 
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3D  RSS  Error  =  32.75  m 


Figure  5.9  -  Simulation  2  RIC  DMC  Results 


Table  5.2  -  Summary  of  Results  with  4  TrakSats  and  50  m  Noise  on  Range 


Simulation  2 

Error  Mean  (m) 

Error  RMS  (m) 

R 

I 

■i,:---.  .<3  ■: 

R 

I 

C 

3D 

No  Process  Noise 

-8.1 

25.8 

-22.4 

53.7 

98.8 

126.5 

169.3 

GEODE  Process  Noise 

-1.4 

7.04 

-2.9 

28.7 

39.5 

12.1 

50.3 

XYZ  DMC 

-2.1 

0.9 

-6.4 

19.1 

26.1 

13.7 

35.2 

RIC  DMC 

-3.5 

2.6 

-4.9 

15.5 

25.9 

12.6 

32.8 

5.6  Simulation  Conclusions 

In  both  simulations  the  RIC  implementation  of  DMC  outperforms  the 
XYZ  version.  The  poor  viewing  geometry  in  the  first  simulation  contributes  to 
allowing  RIC  DMC  to  best  XYZ  DMC  by  over  36%.  The  second  simulation 
shows  that  with  better  geometry  and  more  observations  (even  though  50  m  noise 
is  added)  RIC  DMC  does  not  outperform  XYZ  DMC  as  dramatically.  The 


improvement  in  the  second  simulation  is  slightly  under  7%.  The  several  meter  3D 
improvement  is  worth  the  difficulty  of  tuning  6  parameters  rather  than  2.  The 
greater  flexibility  offered  by  RIC  DMC  in  adjusting  the  process  noise  allowed  a 
19%  improvement  in  the  radial  direction.  Next  results  from  DMC 
implementation  in  GEODE  are  presented. 

5.7  XYZ  DMC  in  GEODE 

The  XYZ  implementation  of  DMC  in  GEODE  is  straightforward.  The 
three  compensative  accelerations  are  added  to  the  state  and  the  deterministic 
portion  of  equation  (5.6)  (also  shown  in  equation  (5.23))  augments  GEODE’s 
time  propagation  of  the  state.  The  state  transition  matrix  is  augmented  with  three 
rows  and  three  columns  per  equation  (5.24).  The  process  noise  covariance 
formulation  in  equations  (5.36)  through  (5.43)  is  implemented  in  the  same  manner 
as  GEODE’ s  process  noise  formulation.  Finally,  x,  ctu  ,  cxb ,  cr6 ,  ctCd  ,  and  ctCr  input 

variables  are  added  to  the  uplink  command  file  along  with  a  switch  to  use  DMC 
or  GEODE’s  original  process  noise  implementation. 

5.7.1  XYZ  DMC  Tuning 

Significant  effort  is  involved  in  tuning  the  process  noise  x  and  o  s  in  the 
application  of  DMC.  The  tuning  process  consists  of  incrementally  changing  each 
parameter  and  recording  the  error  between  the  filter’s  estimate  and  the  “truth” 
data.  An  example  of  a  series  of  tuning  runs  for  the  4  Feb  1997  GPS/MET  data  is 


presented  in  Table  5.3.  In  Table  5.3  the  bold  items  are  the  values  being  tuned. 
The  parameters  that  achieved  optimal  results  are  presented  in  Table  5.4.  The 


results  achieved  with  XYZ  DMC  are  a  significant  improvement  over  GEODE’ 
original  process  noise  formulation. 


Table  5.3  -  GPS/MET  XYZ  DMC  Parameter  Tuning 


T 

<*u 

Gb 

ar 

ac 

'"R 

3DRSS 

0.028 

0.0001 

11.0 

0.001 

0.0001 

0.0001 

8.33 

0.03 

0.0001 

11.0 

0.001 

0.0001 

0.0001 

8.32 

0.034 

0.0001 

11.0 

0.001 

0.0001 

0.0001 

8.33 

0.03 

0.00005 

11.0 

0.001 

0.0001 

0.0001 

8.49 

0.03 

0.00013 

11.0 

0.001 

0.0001 

0.0001 

8.35 

0.03 

0.0001 

10.0 

0.001 

0.0001 

0.0001 

8.34 

0.03 

0.0001 

12.0 

0.001 

0.0001 

0.0001 

8.34 

0.03 

0.0001 

11.0 

0.0001 

0.0001 

0.0001 

8.32 

0.03 

0.0001 

11.0 

0.01 

0.0001 

0.0001 

8.36 

0.03 

0.0001 

11.0 

0.001 

0.00001 

0.0001 

8.32 

0.03 

0.0001 

11.0 

0.001 

0.0002 

0.0001 

8.34 

0.03 

0.0001 

11.0 

0.001 

0.0001 

0.00001 

8.32 

0.03 

0.0001 

11.0 

0.001 

0.0001 

0.0002 

8.34 

Table  5.4  -  GPS/MET  XYZ  DMC  “Optimal”  Parameters 


Parameter 

Tuned  Value 

T 

0.03 

0.0001 

tfb 

11.0 

Gb 

0.001 

acD 

0.0001 

Gcr 

0.0001 

Optimal  tuning  results  for  GPS/MET  data  when  JPL  high  rate  GPS  clock 
estimates  are  applied  to  attempt  to  remove  SA  affects  are  shown  in  Table  5.5. 


Table  5.5  -  GPS/MET  With  High  Rate  GPS  Clocks  XYZ  DMC  “Optimal” 

Parameters 


Parameter 

Tuned  Value 

X 

0.03 

0.0001 

<*b 

5.9 

0.09 

acD 

0.0001 

acR 

0.0001 

5.7.2  XYZ  DMC  GEODE  Accuracy/Precision  Results 


XYZ  DMC  implementation  significantly  improves  GEODE’s  position 
estimates.  Table  5.6  shows  a  comparison  of  the  GPS/MET  (no  high  rate  clocks) 
results  achieved  with  GEODE’s  original  process  noise  formulation  and  results 
with  XYZ  DMC  in  GEODE.  Without  DRVID  XYZ  DMC  improves  the  error 
RSS  by  over  23.0%  while  DRVID  XYZ  DMC  improves  accuracy  by  over  28.2%. 
Both  of  the  previous  comparisons  are  with  the  original  implementation  of 
GEODE. 


Table  5.6  -  GEODE  Results  With  and  Without  XYZ  DMC 


4  Feb  1997  GPS/MET 

Mean  of  Error  (m) 

RMS  Error  of  Filter  Compared 
to  Truth  (m) 

Scheme 

R 

I 

c 

R 

I 

c 

3D 

Original  Single 

0.55 

-0.99 

-0.10 

2.74 

10.24 

2.15 

10.81 

Original  DRVID 

0.24 

-1.14 

-0.11 

2.63 

9.72 

2.23 

10.32 

XYZ  DMC  Single 

-0.18 

3.08 

-0.14 

1.74 

7.77 

2.41 

8.32 

XYZ  DMC  DRVID 

-0.30 

1.94 

-0.14 

1.93 

7.07 

2.56 

7.76 

To  gauge  the  accuracy  improvement  gained  using  XYZ  DMC  with  SA  off, 
results  are  presented  in  Table  5.7  and  Table  5.8.  Table  5.7  shows  DMC  results 


when  JPL  high  rate  GPS  clock  estimates  are  applied  to  GPS/MET  data.  Table  5.8 
shows  DMC  results  for  the  S  A  free  TOPEX  data.  DRVID  XYZ  DMC  provided  a 
19.6%  improvement  over  the  original  implementation  in  GEODE  while  XYZ 


DMC  without  DRVID  resulted  in  only  a  16.7%  improvement  on  the  GPS/MET 


data.  Also,  DMC  improves  the  SA  free  TOPEX  3D  position  RSS  by  10.8%  when 


precise  GPS  ephemerides  are  not  used  and  11.1%  when  precise  GPS  ephemerides 


are  used. 


Table  5.7  -  GEODE  Results  With  and  Without  XYZ  DMC  Using  High  Rate 
Clock  Estimates  to  Correct  for  S  A 


4  Feb  1997  GPS/MET 

Mean  of  Error  (m) 

RMS  Error  of  Filter  Compared 
to  Truth  (m) 

Scheme 

R 

I 

R 

■  I 

c 

3D 

Original  None 

0.28 

-3.59 

-0.04 

2.22 

7.23 

0.77 

7.60 

Original  DRVID 

-0.03 

-3.54 

-0.04 

1.97 

6.95 

0.90 

7.28 

XYZ  DMC  None 

-0.21 

-0.52 

-0.05 

2.01 

5.94 

0.89 

6.33 

XYZ  DMC  DRVID 

-0.31 

-1.19 

-0.06 

1.95 

5.70 

1.03 

6.11 

Table  5.8  -  SA  Free  TOPEX  Results  With  and  Without  XYZ  DMC 


5  May  2000  TOPEX 

Mean  of  Error  (m) 

RMS  Error  of  Filter  Compared 
to  Truth  (m) 

Scheme 

R 

I  V 

C  : 

R 

I 

C 

3D 

Original  None 

0.21 

-0.04 

0.01 

0.35 

1.06 

0.44 

1.20 

XYZ  DMC  None 

0.19 

-0.03 

0.02 

0.33 

0.91 

0.45 

1.07 

Original  Single  with 
Precise  GPS 

Ephemeris 

0.22 

-0.16 

0.00 

0.34 

0.97 

0.31 

1.08 

XYZ  DMC  Single 
with  Precise  GPS 
Ephemeris 

0.20 

-0.18 

0.01 

0.33 

0.84 

0.33 

0.96 
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5.7.3  XYZ  DMC  GEODE  Computational  Burden  Results 

The  computational  burden  of  XYZ  DMC  is  roughly  equivalent  to  that  of 
the  original  process  noise  implementation  in  GEODE. 

5.8  RIC  DMC  in  GEODE 

The  implementation  of  RIC  DMC  in  GEODE  involves  the  same 
implementation  as  XYZ  DMC  with  the  addition  of  the  RIC  variables  to  the  uplink 
command  file  and  appropriate  application  of  the  XYZ  to  RIC  transformation 
matrices  as  in  equations  (5.46),  (5.47),  (5.54),  and  (5.71).  It  can  be  shown  that  by 
setting  tR  =  Tj  =  xc  and  aR  =  Cj  =  ac  RIC  DMC  is  the  equivalent  of  XYZ  DMC. 

The  first  line  of  data  in  Table  5.9  shows  the  XYZ  DMC  case  applied  to  RIC 
DMC.  Table  5.9  also  shows  the  results  achieved  with  GEODE  through  tuning  of 
the  RIC  parameters  for  the  case  where  single  frequency  GPS  measurements  are 
used  and  DRVID  is  not  applied.  Table  5.10  shows  the  optimal  RIC  DMC 
parameters  found  through  tuning. 


Table  5.9  -  GPS/MET  RIC  DMC  Parameter  Tuning 


Tc 

3D  RSS 

0.03 

0.03 

0.03 

0.0001 

0.0001 

0.0001 

8.32 

0.001 

0.03 

0.03 

0.0001 

0.0001 

0.0001 

8.32 

0.1 

0.03 

0.03 

0.0001 

0.0001 

0.0001 

8.35 

0.03 

0.028 

0.03 

0.0001 

0.0001 

0.0001 

8.33 

0.03 

0.032 

0.03 

0.0001 

0.0001 

0.0001 

8.33 

0.03 

0.03 

0.3 

0.0001 

0.0001 

0.0001 

8.10 

0.03 

0.03 

0.4 

0.0001 

0.0001 

0.0001 

8.08 

0.03 

0.03 

0.5 

0.0001 

0.0001 

0.0001 

8.08 

0.03 

0.03 

0.6 

0.0001 

0.0001 

0.0001 

8.09 

0.03 

0.03 

0.4 

0.00005 

0.0001 

0.0001 

8.07 

0.03 

0.03 

0.4 

0.0002 

0.0001 

0.0001 

8.09 

0.03 

0.03 

0.4 

0.00005 

0.00005 

0.0001 

8.07 

0.03 

0.03 

0.4 

0.00005 

0.0002 

0.0001 

8.09 

0.03 

0.03 

0.4 

0.00005 

0.0001 

0.00005 

8.07 

0.03 

0.03 

0.4 

0.00005 

0.0001 

0.0002 

8.09 

Table  5.10  -  GPS/MET  RIC  DMC  “Optimal”  Parameters 


jy  Parameter 

Tuned  Value 

TR 

0.03 

*1 

0.03 

XC 

0.4 

<*R 

0.00005 

CTI 

0.0001 

CTC 

0.0001 

°b 

11.0 

CTb 

0.001 

acD 

0.0001 

acR 

0.0001 

5.8.1  RIC  DMC  GEODE  Results 

Table  5.1 1  shows  a  comparison  of  the  RIC  DMC  results  using  the  optimal 
values  shown  in  Table  5.10  compared  to  the  original  process  noise 
implementation  in  GEODE.  The  position  error  improved  from  8.32  m  3D  RSS 


for  XYZ  DMC  to  8.07  m  3D  RSS  for  RIC  DMC  or  3%  (case  without  DRVID). 
The  DRVID  improvement  is  only  2.6%. 

The  accuracy  improvements  realized  by  RIC  DMC  in  GEODE  are  much 
smaller  than  those  gained  by  RIC  DMC  in  the  simulations  previously  described. 
The  small  improvement  in  accuracy  is  probably  due  to  the  exceptional  observing 
geometry  offered  by  GPS.  As  shown  in  the  simulations,  when  the  observing 
geometry  improves,  the  utility  of  RIC  DMC  decreases.  Although  RIC  DMC 
implementation  only  improves  accuracy  marginally  for  the  processing  of  GPS 
pseudoranges,  it  may  provide  additional  accuracy  improvements  when  applied  to 
different  observation  types,  as  in  the  case  where  GPS  point  solutions  are  used  as 
measurements  for  a  filter.  Additionally,  RIC  DMC  increased  GEODE’s 


computational  burden  by  approximately  1%. 


Table  5.11  -  GEODE  Results  With  and  Without  RIC  DMC 


4  Feb  1997  GPS/MET 

Mean  of  Error  (m) 

RMS  Error  of  Filter  Compared 
to  Truth  (m) 

Scheme 

R 

I 

C 

R 

V  I 

C 

3D 

Original  Single 

0.55 

-0.99 

-0.10 

2.74 

10.24 

2.15 

10.81 

Original  DRVID 

0.24 

-1.14 

-0.11 

2.63 

9.72 

2.23 

10.32 

XYZ  DMC  Single 

-0.18 

3.08 

-0.14 

1.74 

7.77 

2.41 

8.32 

XYZ  DMC  DRVID 

-0.30 

1.94 

-0.14 

1.93 

7.07 

2.56 

7.76 

RIC  DMC  Single 

-0.17 

2.89 

-0.17 

1.70 

7.55 

2.30 

8.07 

RIC  DMC  DRVID 

-0.29 

1.75 

-0.17 

1.94 

6.92 

2.26 

7.53 

5.9  Summary 


In  this  chapter  Dynamic  Model  Compensation  is  first  developed  for  a  one¬ 
dimensional  estimation  problem.  DMC  is  then  extended  to  a  satellite  orbit 


determination  problem  and  the  XYZ  and  RIC  versions  are  developed.  XYZ  and 
RIC  DMC  are  then  employed  in  two  simulations.  In  the  first  simulation,  RIC 
DMC  lowers  the  3D  error  RSS  by  36.4%  compared  to  XYZ  DMC.  Both  XYZ 
and  RIC  DMC  improved  accuracy  considerably  over  GEODE’s  process  noise 
implementation.  In  the  second  simulation,  observing  geometry  is  improved 
causing  the  RIC  DMC  accuracy  to  only  improve  7%  over  XYZ  DMC. 

XYZ  and  RIC  DMC  are  then  implemented  in  GEODE.  XYZ  DMC 
improves  GEODE’s  error  RSS  by  23.0%  when  applied  to  the  processing 
GPS/MET  data,  16.7%  for  GPS/MET  data  with  high  rate  GPS  clocks  and  10.8% 
on  SA  free  TOPEX  data.  XYZ  DMC  not  only  improves  the  position  error  RSS,  it 
does  not  increase  computational  burden. 

While  RIC  DMC  significantly  improves  filter  accuracy/precision  in  the 
simulations  it  only  marginally  improves  GEODE’s  performance.  The  modest  3% 
improvement  without  DRVID  and  2.6%  improvement  with  DRVID  (compared  to 
XYZ  DMC)  may  not  be  worth  the  overhead  involved  in  tuning  the  additional 
parameters.  Finally,  RIC  DMC  increases  GEODE’s  computational  burden  by 
roughly  1%  compared  to  GEODE’s  original  process  noise  formulation. 


CHAPTER  6 


Genetic  Model  Compensation  (GMC) 

Genetic  Model  Compensation  (GMC)  [9,  64]  is  an  application  of  a 
Genetic  Algorithm  (GA)  to  the  optimization  of  the  correlation  time  of  the 
exponential  decay  (x)  and  standard  deviation  (au,  ab,  at )  of  the  forcing  noise 

from  DMC.  As  seen  in  Chapter  5  of  this  dissertation,  the  values  chosen  for  x  and 
the  as  have  a  significant  effect  on  filter  performance.  Normally,  the  as  are 
considered  constant  and  tuned  through  trial  and  error  testing  where  the  true  states 
are  available  (as  shown  in  Chapter  5).  Since  x  is  included  in  the  deterministic 
portion  of  the  DMC  augmentation  of  the  state,  it  can  be  estimated  as  part  of  the 
state.  Unfortunately,  there  are  problems  associated  with  estimating  x.  The 
satellite  OD  EKF  relies  on  linearization  of  the  state  trajectory  and 
observation/state  relationship.  The  state  members  (position,  velocity,  clock  bias, 
clock  bias  drift  rate,  coefficient  of  drag  and  coefficient  of  radiation  pressure)  are 
known,  a  priori,  with  sufficient  accuracy  to  guarantee  linearity.  However,  there  is 
no  way  to  generate  an  a  priori  estimate  of  x  with  accuracy  to  support  the 
assumption  of  linearity.  Without  accurate  a  priori  knowledge  of  x,  the  filter  may 
not  be  able  to  adjust  x  to  its  optimum  value  and  may  even  diverge  [65].  If  x  and 
the  as  are  tuned,  there  is  no  guarantee  their  optimum  value  will  be  chosen  or  the 


optimum  tuned  values  will  be  optimum  during  actual  implementation  (flight  on 
board  a  satellite,  in  this  case).  GMC  was  developed  and  implemented  by 
Cruickshank  to  replace  ad  hoc  tuning  and  estimation  with  the  genetic  estimation 
of  au,  ab,  crb  and  x  [9, 64, 65].  GMC  adapts  au,  ab,  a6  and  x  as  part  of  the  OD 

process.  Cruickshank  suggests  that  GMC  is  perfectly  suited  for  the  autonomous 
satellite  OD  environment  since  no  a  priori  knowledge  of  ctu  ,  ab ,  ab  and  x  is 

needed  [9]. 

In  this  chapter  the  Genetic  Algorithm  (GA)  is  briefly  described,  the 
algorithm  for  the  application  of  GMC  in  GEODE  is  presented  and 
accuracy/precision  and  computational  burden  results  of  GMC  are  shown. 

6.1  The  Genetic  Algorithm  (GA) 

A  summary  of  the  benefits  of  the  GA  and  a  brief  description  are  presented 
here,  more  detailed  information  can  be  found  in  [66].  The  GA  is  a  model  of 
machine  learning  derived  from  the  theoretical  mechanisms  of  evolution  in  nature. 
In  the  GA  chromosomes  represent  a  population  of  individuals.  The  chromosomes 
are  a  set  of  binary  numbers  representing  values  to  be  optimized.  The  individuals 
go  through  a  process  of  simulated  evolution  to  optimize  their  values  based  on  a 
fitness  function.  There  are  many  optimization  techniques  available  however,  the 
GA  has  several  qualities  that  separate  it  from  other  search  techniques.  The  GA: 

•  Operates  on  coding  of  the  parameter  set  rather  than  on  the  parameters 
themselves 


•  Searches  for  the  optimum  from  a  population  of  points  simultaneously 
rather  than  from  a  single  point 

•  Uses  payoff  information  from  the  values  generated  by  a  fitness  function 
rather  than  the  gradient  or  other  auxiliary  information  to  direct  the  search 

•  Uses  probabilistic  transition  rules  rather  than  deterministic  transition  rules 

[66] 

The  GA  seeks  a  global  maximum  (or  minimum)  by  optimally  using 
exploration  to  investigate  unknown  areas  of  the  search  space  and  exploitation  to 
make  use  of  knowledge  found  at  previously  visited  points  [67],  The  three 
operators  reproduction,  crossover  and  mutation  act  on  the  population  allowing 
both  exploration  and  exploitation  of  the  search  space. 

6.1.1  Reproduction 

In  reproduction,  parents  are  selected  randomly  from  the  population  in  a 
way  that  favors  the  most  sought  after  individuals.  The  most  sought  after 
individuals  may  be  selected  more  than  once,  while  least  sought  after  individuals 
may  not  be  selected.  The  parents  (chosen  for  their  fitness  for  reproduction  based 
on  a  fitness  function)  are  then  randomly  arranged  in  pairs  for  mating. 

6.1.2  Crossover 

Crossover  is  the  random  exchange  of  chromosome  strings  between  mating 
pairs.  Crossover  is  not  always  performed  on  the  entire  population  but  on  random 
pairs  of  individuals  selected  for  mating.  The  probability  of  selection  for  crossover 
is  between  50  and  100%.  Random  portions  of  the  chromosomes  from  one 


individual  are  exchanged  for  those  from  another.  In  this  way  the  offspring  may 
receive  a  selection  of  genes  from  each  parent. 

6.1.3  Mutation 

Mutation  is  randomly  applied  to  offspring  after  crossover.  Here  individual 
genes  are  randomly  altered  with  very  small  probability  (0.1  -  0.001).  The 
traditional  view  is  that  crossover  is  more  important  for  rapidly  exploring  the 
search  space  while  mutation  provides  a  small  amount  of  random  search,  ensuring 
that  no  point  in  the  search  space  has  zero  probability  of  examination  [67]. 

6.1.4  GA  Application 

The  GA  is  applied  by,  first,  randomly  selecting  a  population  of 
individuals.  These  individuals  are  tested  by  a  fitness  function  to  determine  which 
are  most  suitable  for  reproduction.  Then,  a  probability  of  reproduction  is 
calculated  for  individuals  based  on  their  suitability.  Reproduction  is 
accomplished  by  randomly  selecting  individuals  based  on  their  reproduction 
probability.  In  this  way,  the  best  individuals  have  the  highest  priority  of  being 
selected.  Individuals  may  be  selected  more  than  once  or  not  at  all.  Next, 
individuals  are  randomly  paired  for  mating.  Then,  crossover  of  genes  is 
performed  if  randomly  decided.  Finally,  mutation  is  performed  on  offspring  from 
crossover,  again  if  randomly  decided  [64]. 
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6.2  GMC  Algorithm 

The  following  GMC  algorithm  is  extracted  from  Cruickshank  [9]  and 
applies  to  the  pseudorange  measurement  only.  The  algorithm  is  modified  to 

include  <x  . 

0 

1.  Generate  an  initial  gene  pool  (as  base  10  numbers)  of  size  n  for 

au ,  ob ,  ab  and  x.  Nominally,  n  =10.  A  sample  distribution  for  the 

gene  pool  is  [1,  50, 100,  500, 1000,  3300, 6600, 10000,  13000,  16383]. 
Scaling  of  the  gene  pool  is  discussed  in  section  6.3. 

2.  Encode  the  n  base  ten  numbers  (candidate  values  for  <xu,  ab,  crb,  x)  in 

each  gene  pool  into  n  binary  strings  of  length  q  (nominally,  q=14). 

3.  Randomly  group  the  candidate  values  into  quadruplets  [x  <ru  ab  ab  J . 

4.  Compute  the  propagated  genetic  position,  velocity,  and  clock  increments 

for  each  of  the  n  candidate  quadruplets  [x  au  ab  ab  ]  using: 


5xGMC  =  f‘  [x(t-T) + x2  (e-(t'T)/T  -l)u(t)dT]  (6.1) 

8b=  f  '  ( t-T  )ub  (T)dT  (6.2) 

5bu=fub(T)dT  (6.3) 

Jt0 

The  stochastic  integrals  are  evaluated  by  Riemann  sums,  i.e.: 

m 

f  g(u(T))dTs £g(u(t;))(tw-t,)  (6.4) 

0  i=0 

m  (  At  ^ 

5xGmc=Z  tksx  +  x2  (eks  - 1) — u  (t)ou  (6.5) 

j=i  V  m  J 

m  (  A f  At  ^ 

5b  =]£  tks— ub  (t)ob  + — ub(i)as  (6.6) 

hI  m  m  J  J 


At 

—  Ubj(t)CTb 
m  1 
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8b  =£ 

j=i 


tks  =  At 


(m- 0.5)  At 


eks  =  e 


where  t;  <  t;  <  ti+1  and  u(t)  represents  the  uncorrelated  Gaussian  random 
noise  process. 

5.  Compute  the  predicted  genetic  pseudorange  by 

°r!  =  Pi  -  b.c  +  b<=c  =  J(xG-x«)2-(yc-y,)2-(zG-z,)2  -  b,c  +  bGc  (6. 10) 


where 


xg=xu+5xGMC 


yg  =  yu+5yc 


(6.11) 


zg  —  zu  +  8zgmc 


bg  —  bu  +  5bGMC  +  bu  At  +  8bGMCAt 


bo  =  GPS  satellite  clock  bias 

xu,  yu,  and  zu  are  the  best  estimates  of  the  user  satellite  at  the  current  epoch 
and  xg,  yG,  and  zG  are  the  GPS  satellite’s  coordinates. 


6.  Compute  the  genetic  pseudorange  residuals  by 

y  PEi  —  Pi  —  ^pgi 

where  p*  is  the  ith  pseudorange  measurement. 

7.  Compute  the  fitness  function  value  for  each  candidate  quadruplet 

[x  au  ob  ab]  by 


(6.12) 


fi=7T^T-;i  =  1’2’  •n 

£g(K  8gi 


(6.13) 


where 
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i  =  the  ith  genetic  parameter  quadruplet 

sgi  =  the  vector  of  genetic  measurement  residual  for  the  ith  genetic 
parameter  quadruplet 

R'1  =  a  diagonal  weighting  matrix  formed  from  the  observed  measurement 
residuals 


a  =  the  maximum  value  of  sg  R  le  over  all  n  genetic  parameter 
quadruplets 

8.  Select  the  candidate  quadruplet  [x  <tu  ab  ct^  ]  that  yields  the 
highest  objective  function  value: 

^max  ~  f  (Topt » auop(  ’  abop,  >  abopt  ) 


(6.14) 


9.  Update  |^x  cru  ab  ab  J  (  for  the  next  formulation  of  the  process  noise 

covariance  and  state  transition  matrices. 

10.  Update  the  propagation  equations  for  the  acceleration,  velocity,  and 

position  states  so  that  Topt  will  be  used  in  the  next  time  update. 

11.  Compute  the  merit  values  (probability  of  reproduction)  for  the  gene  pool 

quadruplets: 

v,  =TT—  (6.15) 

If, 

j=l 


12.  Compute  the  gene  pool  variances.  For  any  particular  gene  pool,  the  gene 
pool  variance  is  defined  as: 

«8)  =  1Z(9,-ew)!  (6.16) 

n  i=l 


where  0  represents  the  gene  pool  parameter,  e.g.,  a  or  x. 

13.  Perform  the  reproduction  procedure.  The  probability  of  reproduction  for 
the  ith  gene  pool  quadruplet  is  given  by  Vj.  Reproduction  of 
au,  Qb,  ab  and  x  are  done  by  quadruplets. 


14.  Perform  the  crossover  procedure  with  the  new  gene  pools.  The  crossover 
probability  is  a  function  of  the  gene  pool  variance: 
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¥(0)  =  0.99  -  0.99 


K(0) 


Pcross 

Pcross 


(0)  =  VF(0)  if  ^(0)>O.5 
(0)  =  0.5  if  'F(0)<O.5 


(6.17) 


where  k(0)  is  the  crossover  rate  factor.  Its  value  is  chosen  to  be  of  the 
same  order  of  magnitude  as  the  anticipated  maximum  value  of  ^(0) . 
The  probability  floor  of  0.5  is  arbitrarily  chosen.  Crossover  is  performed 
on  the  au ,  CTb ,  ,  t  strings  separately,  not  on  the  quadruplets. 


15.  Perform  the  mutation  procedure  with  the  new  gene  pools.  The  probability 

of  mutation  for  any  particular  bit  is  an  input  value  (nominally  set  to 
0.01).  Mutation  is  performed  on  the  ctu,  ab,  x  strings  separately. 

16.  Return  to  step  3  with  the  new  au,  crb,  o6,  x  gene  pools.  The  new  gene 

pools  become  the  initial  gene  pools  for  the  next  update  at  tk+i- 


Figure  6.1  shows  the  entire  GMC  algorithm  with  respect  to  GEODE. 
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Figure  6.1  -  GMC  Algorithm 


6.3  GMC  in  GEODE 

GMC  is  implemented  in  GEODE  using  XYZ  DMC  only.  The  values  for 
<jcd  and  ocr  are  set  at  0.0001,  just  as  in  XYZ  DMC.  The  binary  string  length  is 

set  at  14  bits  giving  a  range  of  values  for  each  member  of  the  design  space  from 
0  to  16383  (214).  Each  x  and  ct  is  scaled  to  allow  the  design  space  to  more 
appropriately  match  the  tuned  DMC  values.  Additional  information  concerning 
the  constraints  placed  on  the  design  space  for  each  t  and  a  is  presented  in  section 
6.3.2.  The  size  of  the  population  of  quadruplets  is  set  to  10  (n  =  10);  therefore. 


there  are  40  members  in  the  design  space.  Fifty  intervals  are  used  for  the 
Riemann  sum  calculations  shown  in  equations  (6.5)  through  (6.7),  i.e.,  m  =  50. 
The  population  size  and  number  of  Riemann  intervals  are  changeable  parameters. 

6.3.1  Implementation 

The  implementation  of  GMC  in  GEODE  is  very  straightforward.  The 
previous  implementation  of  XYZ  DMC  is  left  untouched  and  the  part  of  Figure 

6.1  inside  the  dashed  lines  is  added  to  GEODE.  All  of  the  GMC  processes  take 
place  just  after  the  time  update  of  the  state  and  just  prior  to  the  formulation  of  the 
state  transition  and  process  noise  covariance  matrices. 

6.3.2  GMC  Tuning 

Although  trial  and  error  tuning  is  not  supposed  to  be  required  in  GMC  [9], 
GMC  results  improve  dramatically  if  the  range  of  values  allowed  for  each 
member  of  the  design  space  is  tightly  constrained  around  the  tuned  DMC  values. 
This  “bracketing”  of  the  design  space  is  carried  out  when  the  binary  strings  are 
converted  to  real  numbers.  The  values  chosen  to  bracket  the  design  space  for 
both  the  GPS/MET  data  with  high  rate  GPS  clocks  applied  and  SA  free  TOPEX 
data  are  shown  in  Table  6.1.  Table  6.1  also  shows  the  best-tuned  DMC  values. 

Scaling  of  the  Riemann  sums  is  also  performed  to  keep  the  genetic 
measurement  residuals  within  reasonable  values.  Scaling  of  the  Riemann  sums  is 
performed  in  a  very  straightforward  and  logical  manner.  Once  accomplished,  the 


Riemann  sum  scaling  values  are  used  for  all  data  sets  processed  with  GEODE. 


The  8x,  5y,  and  8z  Riemann  sums  are  scaled  (multiplied)  by  le6  and  the  8b  and 


the  6bdot  sums  are  scaled  by  100. 


Table  6.1-  GMC  Design  Space  Bracketing 


Data  Set 

Best  DMC  Value 

Lower  Bound 

Upper  Bound 

4  Feb  1997  GPS/MET 
with  High  Rate  GPS 
Clock  Estimates 

t 

0.03 

0.001 

0.05 

au 

0.0001 

0.00005 

0.00015 

ab 

5.9 

0.001 

8.0 

tfbdot 

0.09 

0.01 

0.4 

5  May  2000  TOPEX 

T 

0.01 

0.001 

0.05 

CTu 

0.0001 

0.00005 

0.00015 

0.01 

0.001 

8.0 

^bdot 

0.2 

0.01 

0.4 

6.3.3  GMC  in  GEODE  Results 

The  application  of  the  GA  in  GMC  necessitates  the  use  of  a  significant 
quantity  of  random  numbers.  GMC  relies  on  random  numbers  for  decisions  in 
reproduction,  crossover  and  mutation  as  well  as  in  the  solution  of  the  stochastic 
integrals  using  Riemann  sums.  Since  the  generation  of  random  numbers  depends 
on  an  initial  “seed,”  the  results  attained  with  GMC  change  whenever  a  different 
random  number  seed  is  used.  Therefore,  to  gauge  GMC’s  performance  against 
DMC,  statistical  hypothesis  testing  is  implemented.  Hypothesis  testing  is  a 
procedure  where  the  mean  of  a  population  can  be  compared  against  a  specified 
value  [53].  In  this  case,  the  position  RIC  RMS  and  RSS  values  achieved  using 


DMC  and  JPL  high  rate  GPS  clock  estimates  (GPS/MET  data)  are  compared  to 
the  results  of  30  similar  GMC  runs.  The  sample  size  of  30  runs  is  chosen  to 
satisfy  the  Central  Limit  Theorem.  The  Central  Limit  Theorem  states  that  if  the 
sample  size  n  is  large  (n  >  30) ,  the  sampling  distribution  of  the  sample  mean  will 
be  approximately  normal,  even  if  the  probability  distribution  of  the  population  is 
unknown  [53].  GEODE  is  run  using  GMC  30  times  with  random  random  number 
seeds.  The  radial,  in-track,  cross-track  position  errors  from  the  30  runs  are  used 
in  the  hypothesis  testing. 

The  hypothesis  testing  procedure  as  shown  in  [53]  is  set  up  below  to 
compare  GMC  applied  to  the  GPS/MET  data  with  high  rate  GPS  clock  estimates 
against  GEODE’s  original  process  noise  formulation  with  high  rate  GPS  clock 
estimates: 

1.  The  parameter  of  interest  is  the  mean  of  the  position  RSS  of  the  30  GMC 
runs,  referred  to  as  p. 

2.  The  null  hypothesis  is  Ho:  p  =  7.60  m  (3D  position  RSS  of  the  original 
process  noise  implementation  in  GEODE) 

3.  The  alternative  hypothesis  is  Hi :  p  <  7.60m  (reject  Ho  if  the  mean  of  the 
position  RSS  of  the  GMC  runs  is  less  than  7.60  m) 

4.  The  type  I  error  probability  (the  probability  of  rejecting  the  null 
hypothesis  even  when  it  is  true)  is  set  at  a  =  0.0001. 

5.  The  test  statistic  is  t0  =  — — (6. 

s/>/n 

6.  Reject  Ho  if  to  <  to.oooi,  29  =  4.254 

7.  Computations:  Since  x  =  6.13  m  (sample  mean),  s  =  0.317  (sample 

standard  deviation),  po  =  7.60  m  and  n  =  30, 
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t0 


6.13-7.6 

0.317/^0 


-28.79 


(6.19) 


8.  Conclusions:  Since  to  =  -28.79  <  4.254,  we  reject  Ho  and  conclude,  at  the 

0.0001  significance  level,  the  mean  of  the  GMC  position  RSS  is  lower 
than  the  DMC  position  RSS.  In  other  words,  there  is  a  99.99%  probability 
GMC  produces  a  lower  3D  position  RSS  than  DMC. 


Similar  hypothesis  tests  are  performed  for  the  radial,  in-track  and  cross¬ 
track  RMS  values  with  the  result  that,  at  the  0.0001  significance  level,  GMC 
produces  lower  radial  and  in-track  errors  than  GEODE’s  original  implementation. 
The  original  implementation  performs  better  than  GMC  in  the  cross-track 
direction.  Figure  6.2  shows  a  plot  of  mean  GMC  results  and  the  original  process 
noise  implementation  in  GEODE  with  the  JPL  high  rate  GPS  clock  estimates 
applied.  As  reflected  in  the  hypothesis  testing  mentioned  previously,  the  plots  in 
Figure  6.2  show  that  GMC  significantly  outperforms  the  original  process  noise 
implementation  in  GEODE  in  the  radial  and  in-track  directions.  GMC  also 
provides  a  less  biased  estimate  in  the  in-track  direction.  Table  6.2  shows  statistics 
of  the  results  using  the  original  version  of  GEODE,  XYZ  DMC  in  GEODE  and 
GMC  in  GEODE  when  the  JPL  high  rate  GPS  clock  estimates  are  applied. 


GMC  3D  Error  RSS  6.09  m  DMC  3D  Error  RSS  7.60  m 


0  4  8  12  16  20  24 


0  4  8  12  16  20  24 


Figure  6.2  -  GPS/MET  High  Rate  GPS  Clock  GMC/DMC  Comparison 


Table  6.2  -  GEODE  GMC  Results  Comparison  Using  High  Rate  Clock  Estimates 

to  Correct  for  SA 


4  Feb  1997  GPS/MET 

Mean  of  Error  (m) 

RMS  Error  of  Filter  Compared 
to  Truth  (m) 

Scheme 

R 

I 

C 

R 

I 

C 

3D 

Original  single 

0.28 

-3.59 

-0.04 

2.22 

7.23 

0.77 

7.60 

XYZ  DMC  single 

-0.21 

-0.52 

-0.05 

2.01 

5.94 

0.89 

6.33 

XYZ  GMC  single* 

-0.28 

0.42 

-0.6 

1.92 

5.56 

0.81 

6.13 

*  Mean  values  for  30  runs 

The  mean  x  and  o  values  chosen  by  GMC  over  the  30  runs  are  different  than  the 
optimal  values  found  with  DMC.  Figure  6.3  shows  plots  of  the  mean  DMC 
parameters  determined  with  GMC  over  the  30  runs.  The  horizontal  lines  are  the 
bracket  values  shown  in  Table  6.1.  The  upper  bracket  value  for  is  off  the 

plot.  The  mean  values  for  x,  g„,  and  are  close  to  the  tuned  DMC  values.  The 

mean  value  for  ab  is  quite  different  and  indicates  that  additional 
accuracy/precision  gains  may  be  attained  with  further  DMC  tuning. 
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Figure  6.3  -  GPS/MET  Mean  DMC  Parameters  from  GMC 


GEODE  with  GMC  is  also  used  to  process  the  SA  free  TOPEX  data.  The 
mean  results  from  the  thirty  GMC  TOPEX  runs  are  shown  in  Table  5.8  for 
comparison  with  the  original  version  of  GEODE  and  XYZ  DMC  in  GEODE. 
Hypothesis  testing  results  indicate  that  with  99.99%  certainty,  the  means  of  the 
GMC  RSS  results  are  smaller  than  original  GEODE.  Hypothesis  tests  also 
indicate  that  GMC  improves  the  radial  and  in-track  position  estimates,  with 
99.99%  certainty.  Again,  the  original  GEODE  implementation  performed  nearly 
the  same  as  GMC  in  the  cross-track  direction.  Figure  6.4  shows  plots  of  the 
original  GEODE  implementation  and  mean  GMC  results  in  the  radial,  in-track 


and  cross-track  directions.  In  the  radial  and  cross-track  directions  the  GMC 


results  are  nearly  exactly  the  same  as  original  GEODE  results. 


Table  6.3  -  SA  Free  TOPEX  GMC  Results 


5  May  2000  TOPEX 

Mean  of  Error  (m) 

RMS  Error  of  Filter  Compared 
to  Truth  (m) 

Scheme 

R 

I 

c 

R 

I 

c 

3D 

Original  Single 

0.21 

-0.04 

0.01 

0.35 

1.06 

0.44 

1.20 

XYZDMC  Single 

0.19 

-0.03 

0.02 

0.33 

0.91 

0.45 

1.07 

XYZ  GMC  Single* 

0.19 

-0.10 

0.02 

0.34 

0.92 

0.45 

1.08 

*  Mean  values  for  30  runs 
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GMC  3D  Error  RSS  1 .1 3  m  DMC  3D  Error  RSS  1 .07  m 


Figure  6.4  -  TOPEX  GMC/DMC  Comparison 


Figure  6.5  shows  plots  of  the  mean  DMC  parameters  produced  by  GMC. 
In  this  case  only  the  mean  value  of  x  is  significantly  different  from  the  best-tuned 
DMC  value.  Another  interesting  aspect  of  the  plot  is  the  spikes  just  after  hours  8 
and  16.  The  cause  of  the  spikes  is  not  known. 


x  mean  0.002438 


Time  (hours) 


Time  (hours) 


Figure  6.5  -  TOPEX  Mean  DMC  Parameters  from  GMC 


The  computational  cost  of  using  GMC  to  estimate  the  XYZ  DMC 
parameters  is  significant.  Run  time  increases  by  86%  when  GMC  is  used  instead 
of  XYZ  DMC.  While  the  86%  increase  is  significant,  the  overall  run  time  when 
processing  24  hours  of  data  for  GEODE  with  GMC  is  only  90-seconds  on  a  450 
MHz  Pentium  II  with  128  MB  of  RAM.  Also,  the  90-second  execution  time 
includes  a  significant  amount  of  file  i/o  (reading  in  GPS  observations,  GPS 
ephemerides  and  POE  and  writing  out  state  estimates,  covariance  matrices,  error 
estimates  and  the  Kalman  gain). 
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6.4  Summary 

In  this  chapter  the  Genetic  Algorithm  (GA)  as  applied  to  satellite  orbit 
determination  in  the  form  of  Genetic  Model  Compensation  (GMC)  is  presented  as 
an  alternative  to  the  DMC  tuning  process.  While  GMC  does  not  necessarily 
improve  GEODE’s  accuracy/precision  compared  to  DMC,  it  does  provide 
improvement  compared  to  the  original  process  noise  implementation  in  GEODE. 
GMC  also  reduces  the  burden  associated  with  DMC  tuning.  GMC  does  introduce 
a  significant  computational  burden,  however,  it  is  not  so  significant  that  it 
prohibits  GMC  from  potential  deployment  in  a  real-time  satellite  OD  scenario. 
GMC  is  applied  here  with  XYZ  DMC  but  GMC  could  be  extended  to  estimate 
RIC  parameters  associated  with  RIC  DMC. 


CHAPTER  7 


Summary,  Contributions,  Future  Work, 
and  Conclusion 


7.1  Summary 

The  differences  between  post-processing  and  real-time  satellite  orbit 
determination  systems  are  presented  in  Chapter  1.  Details  regarding  JPL’s  GOA- 

n,  NRL’s  OCEAN  and  VMSI’s  MicroCosm®  software  suites  are  provided  as 

examples  of  post-processing  systems.  Real-time  GIPSY  and  an  unnamed  system 
from  the  University  of  Nottingham  are  described  as  real-time  examples.  The 
main  differences  between  post-processing  and  real-time  systems  (with  SA  off)  as 
stated  in  Chapter  1  are: 

•  Post-processing  systems  can  use  the  best/largest  gravity  models  available 
while  real-time  systems  must  use  truncated  models 

•  Post-processing  systems  estimate  GPS  ephemerides  along  with  the 
satellite  state  or  use  precise  GPS  ephemerides  while  real-time  systems  use 
broadcast  GPS  ephemerides 

•  Post-processing  systems  use  sophisticated  ionospheric  estimation  models 
while  real-time  systems  typically  ignore  ionospheric  errors 

•  Post-processing  systems  use  batch  or  batch  sequential  processors  while 
real-time  systems  use  EKFs 


Chapter  2  includes  a  description  of  the  OrbView-1  spacecraft  and  its 
secondary  payload  GPS/MET.  Chapter  2  also  provides  a  detailed  description  of 


GEODE,  details  concerning  the  optimization  of  GEODE,  a  Lagrangian 
interpolation  scheme  applied  to  precise  GPS  ephemerides  and  a  method  for 
partially  removing  the  effects  of  SA  on  the  GPS/MET  GPS  observations. 

Baseline  values  of  GEODE’s  performance  when  processing  GPS/MET  and  T/P 
data  are  also  established  in  Chapter  2.  See  Table  7.1. 

Chapter  3  outlines  the  effect  on  accuracy  of  changing  the  degree  and  order 
of  the  gravity  model  used  in  propagating  satellite  ephemeris  and  in  real-time 
satellite  OD  (see  Table  7.1).  It  is  important  to  again  note  that  GEODE  utilizes  an 
EKF.  Therefore,  the  OD  results  presented  in  Chapter  3  cannot  be  used  to  predict 
the  effect  of  using  truncated  gravity  models  in  a  batch  processing  scenario. 

At  the  GPS/MET  altitude  and  inclination,  there  is  a  4.5  m  accuracy  gain 
when  a  30x30  truncation  of  the  JGM-2  gravity  model  is  used  for  a  24-hour 
propagation  instead  of  a  20x20  truncation  (both  compared  to  a  70x70  model). 
However,  there  is  only  a  1  m  3D  RSS  improvement  in  GEODE’s  OD  estimate 
when  a  30x30  truncation  is  used  as  compared  to  a  20x20  truncation.  While  the 
OD  accuracy  difference  is  small,  the  computation  time  required  for  processing 
one  day  of  GPS/MET  observations  is  1 1%  less  for  the  20x20  truncation  compared 
to  the  30x30  truncation.  Chapter  3  also  shows  there  is  no  advantage  to  using 
different  models  (JGM-3  or  EGM-96)  than  JGM-2  and  that  there  is  only  a  minor 
OD  accuracy/precision  improvement  when  truncations  larger  than  30x30  are  used. 

To  reduce  OD  computational  burden  while  maintaining  accuracy/precision 
the  Gravity  Acceleration  Approximation  Function  (GAAF)  is  also  introduced  in 


Chapter  3.  At  the  GPS/MET  altitude  and  inclination  it  is  found  that  a  20  km 
altitude  range  could  be  fit  with  1st  order  polynomials  to  describe  pseudocenter 
coefficients  at  various  latitudes  and  longitudes.  The  20  km  -  1st  order  polynomial 
GAAF  implementation  requires  only  1.88  MB  of  memory  compared  to  the  4.9 
MB  required  by  the  original  implementation  in  [50].  24-hour  ephemeris 
propagation  and  GEODE  OD  results  confirm  that  accuracy/precision  are 
maintained  when  GAAF  is  used  and  that  the  computational  burden  of  GAAF  is 
equivalent  to  that  of  a  5x5  gravity  model  truncation. 

Measurement  errors  due  to  the  ionosphere  are  the  topic  of  Chapter  4.  The 
ionosphere  is  described  and  the  Differenced  Range  Versus  Integrated  Doppler 
(DRVID)  method  of  estimating  ionospheric  measurement  errors  is  presented. 
Measurement  errors  from  the  ionosphere  are  highest  when  satellites  are  in 
sunlight  and  for  GPS/MET  were  sometimes  as  high  as  50  m  in  1997.  The  RMS 
of  the  dual  frequency  ionospheric  corrections  for  the  4  Feb  1997  GPS/MET  data 
is  2.49  m  (2.75  m  in  sunlight/1.84  m  in  darkness).  Application  of  dual  frequency 
ionospheric  corrections  to  the  GPS/MET  data  only  marginally  improved  the  3D 
position  RSS  and  introduced  an  in-track  bias.  The  reason  for  the  marginal  RSS 
improvement  is  probably  due  to  the  dominating  effects  of  SA  errors. 

DRVID  utilizes  the  difference  in  the  way  the  pseudorange  and  phase 
measurements  are  affected  by  the  ionosphere.  Since  the  ionosphere  advances 
phase  and  delays  range  measurements,  a  linear  combination  of  phase  and  range 
measurements  can  remove  ionospheric  errors  to  first  order.  Unfortunately,  in 


using  the  phase  measurement,  an  unknown  integer  ambiguity  is  introduced.  Since 
estimation  of  phase  integer  ambiguities  in  real-time  is  probably  unrealistic  and  the 
level  of  interaction  required  to  use  an  accurate  model  to  predict  the  initial 
measurement  error  due  to  the  ionosphere  is  probably  unacceptable  in  real-time 
OD  systems,  the  zero  bias  DRVID  method  is  proposed.  Zero  bias  DRVED  takes 
advantage  of  the  fact  that  the  measurement  error  due  to  the  ionosphere  is  smallest 
when  the  GPS  satellite  is  at  its  highest  elevation  angle  with  respect  to  the  user 
satellite.  Since  the  minimum  measurement  error  point  can  be  easily  determined  in 
real-time,  the  measurement  error  due  to  the  ionosphere  at  this  point  is  assumed  to 
be  zero.  All  measurements  prior  to  the  maximum  elevation  point,  in  a  given 
tracking  arc,  are  not  corrected  with  DRVID  while  all  measurements  taken  after 
the  maximum  elevation  point  are  corrected  with  DRVID.  DRVID  introduces  a 
12.1%  computational  burden  increase  to  GEODE  but  improves  the  3D  position 
RSS  when  applied  to  the  four  different  GPS/MET  cases  and  to  the  T/P  data  as 
seen  in  Table  7.2. 

Dynamic  Model  Compensation  (DMC)  is  developed  and  applied  in 
Chapter  5.  DMC  provides  process  noise  formulation  and  estimation  of 
accelerations  not  accounted  for  in  GEODE’s  dynamic  model.  DMC  is  applied  to 
two  Matlab  simulations  with  the  result  that  the  3D  position  error  RSS  of  the 
filtered  solution  compared  to  truth  is  significantly  improved  when  XYZ  DMC  is 
applied  and  improved  further  when  RIC  DMC  is  applied.  However,  when  the 
observation  geometry  is  improved  in  the  second  simulation,  RIC  DMC  does  not 


provide  as  significant  an  improvement  as  in  the  first  simulation  (when  observation 
geometry  is  not  as  good). 

XYZ  DMC  also  improves  GEODE’s  position  estimates.  DMC  parameter 
tuning  results  are  presented  and  it  is  shown  that  XYZ  DMC  improves  GEODE’s 
3D  position  error  RSS  by  23.0%  over  GEODE’s  original  process  noise 
formulation  on  GPS/MET  data  when  high  rate  GPS  clock  estimates  are  not 
applied,  16.7%  when  high  rate  GPS  clock  estimates  are  applied  and  10.0%  on  SA 
free  T/P  data.  RIC  DMC  further  improves  GEODE’s  position  estimates  by  3% 
over  XYZ  DMC.  XYZ  DMC  does  not  increase  GEODE’s  computational  burden 
while  RIC  DMC  increases  the  burden  by  roughly  1%. 

Tuning  of  the  DMC  time  correlation  coefficient  ( x )  and  the  standard 
deviation  of  the  forcing  noise  (a )  can  be  tedious  and  time  consuming.  Also,  the 
pre-mission  tuning  process  may  yield  optimum  values  that  may  not  be  optimum 
during  the  mission  when  “real”  data  is  processed.  Genetic  Model  Compensation 
(GMC)  is  described  and  GMC  GEODE  implementation  results  are  presented  in 
Chapter  6.  GMC  is  an  application  of  a  Genetic  Algorithm  (GA)  to  optimize  the  x 
and  a  s  used  by  DMC.  GA’s  operate  on  strings  of  binary  numbers  and  use 
random  processes  to  arrive  at  optimum  values  through  reproduction,  crossover 
and  mutation. 

Since  random  numbers  drive  GMC,  hypothesis  testing  is  explained  and 
applied  to  determining,  with  known  probability,  if  GMC  improves  GEODE’s 
position  estimates  compared  to  GEODE’s  original  process  noise  implementation. 


Thirty  GEODE  runs  are  performed  using  GMC  on  the  GPS/MET  data  with  high 
rate  GPS  clock  estimates  applied.  The  results  of  the  30  GMC  runs  are  compared 
against  results  from  original  GEODE.  Hypothesis  testing  concludes,  with  99.99% 
probability,  that  GMC  improves  GEODE’s  position  estimates.  In  fact,  the  mean 
GMC  3D  position  error  RSS  is  19.3%  lower  than  original  GEODE.  Thirty  GMC 
GEODE  runs  are  also  performed  on  the  SA  free  T/P  data.  Hypothesis  testing 
again  concludes,  with  99.99%  probability,  that  GMC  GEODE  provides  better 
position  estimates  than  original  GEODE.  In  this  case,  the  GMC  GEODE  3D 
position  RSS  is  10.0  %  lower  than  original  GEODE. 

The  improvements  gained  through  the  application  of  GMC  do  not  come 
without  computational  cost.  GMC  increases  GEODE’s  computational  burden  by 
approximately  86%  but  GEODE’s  relatively  small  overall  burden,  even  with 
GMC,  does  not  prohibit  its  use  in  a  real-time  OD  system. 


Table  7.1  -  Chapter  2  through  Chapter  4  Results  Summary 


Error  Mean  (m) 

Error  RMS  (m) 

Chapter  2 

R 

I 

C 

R 

I 

c 

3D 

GPSMET  orig  single 

0.55 

-0.99 

-0.10 

2.74 

10.24 

2.15 

10.81 

GPSMET  orig  dual 

0.11 

2.04 

-0.11 

2.48 

9.68 

2.12 

10.21 

GPSMET  orig  single  PE 

0.58 

-1.05 

-0.09 

2.99 

10.11 

1.97 

10.73 

GPSMET  orig  dual  PE 

0.11 

2.01 

-0.11 

2.48 

9.22 

2.03 

9.76 

GPSMET  orig  single  hr 

0.28 

-3.59 

-0.04 

2.22 

7.23 

0.77 

7.60 

GPSMET  orig  dual  hr 

-0.08 

-1.38 

-0.04 

2.10 

6.14 

0.86 

6.55 

GPSMET  orig  single  hr  PE 

0.30 

-3.56 

-0.03 

2.31 

6.95 

0.74 

7.36 

GPSMET  orig  dual  hr  PE 

-0.06 

-1.35 

-0.03 

1.91 

5.43 

1.04 

5.85 

T/P  orig  single 

0.21 

-0.04 

0.01 

0.35 

1.06 

0.44 

1.20 

T/P  orig  single  PE 

0.22 

-0.16 

0.00 

0.34 

0.97 

0.31 

1.08 

Chapters 

GPSMET  orig  single 

JGM-2  05x05  (OD) 

1.71 

-5.57 

-0.06 

25.42 

53.87 

29.42 

66.43 

JGM-2  10x10 

0.54 

-2.03 

0.62 

7.66 

18.16 

7.15 

20.96 

JGM-2  20x20 

0.52 

-0.97 

0.02 

3.12 

11.05 

2.93 

11.85 

JGM-2  22x22 

0.66 

-1.14 

0.05 

2.90 

10.53 

2.26 

11.15 

JGM-2  24x24 

0.58 

-0.91 

-0.07 

2.77 

10.56 

2.56 

11.21 

JGM-2  25x25 

0.58 

-0.88 

-0.07 

2.78 

10.36 

2.03 

10.92 

JGM-2  26x26 

0.60 

-1.02 

-0.06 

2.69 

10.18 

2.06 

10.73 

JGM-2  27x27 

0.57 

-0.97 

-0.06 

2.72 

10.26 

2.13 

10.82 

JGM-2  28x28 

0.56 

-0.99 

-0.06 

2.74 

10.31 

2.26 

10.91 

JGM-2  30x30 

0.55 

-0.99 

-0.10 

2.74 

10.24 

2.15 

10.81 

JGM-2  40x40 

0.52 

-0.90 

-0.07 

2.71 

10.14 

2.29 

10.75 

JGM-2  50x50 

0.52 

-0.88 

-0.05 

2.71 

10.12 

2.41 

10.75 

JGM-2  70x70 

0.52 

-0.88 

-0.05 

2.70 

10.11 

2.38 

10.74 

GAAF 

0.52 

-0.88 

-0.05 

2.71 

10.11 

2.38 

10.74 

EGM-96  30x30 

0.55 

-0.97 

-0.11 

2.75 

10.25 

2.18 

10.83 

JGM-3  30x30 

0.52 

-0.91 

-0.11 

2.74 

10.24 

1.98 

10.78 

Chapter  4 

GPSMET  orig  DRVID 

0.22 

-0.86 

-0.12 

2.49 

10.00 

2.24 

10.55 

GPSMET  orig  DRVID  PE 

0.24 

-0.92 

-0.11 

2.59 

9.68 

2.08 

10.24 

GPSMET  orig  DRVID  hr 

-0.03 

-3.54 

-0.04 

1.97 

6.95 

0.90 

7.28 

GPSMET  orig  DRVID  hr  PE 

0.00 

-3.57 

-0.04 

1.92 

6.46 

0.81 

6.78 

T/P  DRVID  single 

0.20 

0.01 

0.01 

0.35 

1.02 

0.45 

1.17 

Table  7.2  -  Chapter  5  and  Chapter  6  Results  Summary 


Chapter  5 

Error  Mean  (m) 

Error  RMS  (m) 

XYZDMC 

R 

I 

c 

R  . 

i 

.  C 

3D 

GPSMET  DMC  single 

-0.18 

3.08 

-0.14 

1.74 

7.77 

2.41 

8.32 

GPSMET  DMC  dual 

-0.34 

5.07 

-0.14 

1.84 

8.44 

2.38 

8.96 

GPSMET  DMC  DRVID 

-0.30 

1.94 

-0.14 

1.93 

7.07 

2.56 

7.76 

GPSMET  DMC  single  PE 

-0.16 

3.09 

-0.12 

2.00 

7.51 

2.19 

8.08 

GPSMET  DMC  dual  PE 

-0.33 

5.07 

-0.12 

1.76 

7.92 

2.22 

8.41 

GPSMET  DMC  DRVED  PE 

-0.26 

2.11 

-0.12 

1.92 

6.92 

2.15 

7.50 

GPSMET  DMC  single  hr 

-0.21 

-0.52 

-0.05 

2.01 

5.94 

0.89 

6.33 

GPSMET  DMC  dual  hr 

-0.35 

0.87 

-0.05 

2.11 

5.97 

0.98 

6.41 

GPSMET  DMC  DRVID  hr 

-0.31 

-1.19 

-0.06 

1.95 

5.70 

1.03 

6.11 

GPSMET  DMC  single  hr  PE 

-0.19 

-0.41 

-0.04 

2.00 

5.23 

0.84 

5.66 

GPSMET  DMC  dual  hr  PE 

-0.33 

0.98 

-0.04 

1.84 

4.92 

1.11 

5.37 

GPSMET  DMC  DRVID  hr  PE 

-0.30 

-1.22 

-0.05 

1.83 

4.82 

1.00 

5.25 

T/P  DMC  single 

0.19 

-0.03 

0.02 

0.33 

0.91 

0.45 

1.07 

T/P  DMC  DRVID 

0.19 

-0.20 

0.02 

0.34 

0.89 

0.46 

1.06 

T/P  DMC  single  PE 

0.20 

-0.17 

0.01 

0.33 

0.84 

0.33 

0.96 

T/P  DMC  DRVID  PE 

0.20 

-0.34 

0.01 

0.33 

0.85 

0.35 

0.98 

RICDMC 

GPSMET  RIC  DMC  single 

-0.17 

2.89 

-0.17 

1.70 

7.55 

2.30 

8.07 

GPSMET  RIC  DMC  DRVID 

-0.29 

1.75 

-0.17 

1.94 

6.92 

2.26 

7.53 

Chapter  6 

GPSMET  GMC  single  hr* 

-0.28 

0.42 

-0.6 

1.92 

5.56 

0.81 

6.13 

T/P  GMC* 

0.19 

-0.10 

0.02 

0.34 

0.92 

0.45 

1.08 

*  Mean  value  for  30  runs 


Results  for  all  GPS/MET  and  T/P  GEODE  runs,  as  described  in  Chapter  2 


through  Chapter  6,  are  presented  in  Table  7.1  and  Table  7.2.  The  following 


definitions  also  apply  to  Table  7.1  and  Table  7.2: 


orig  =  Original  GEODE  process  noise  formulation 
single  s  No  ionospheric  correction  applied 
PE  3  Precise  GPS  ephemerides  used 
hr  s  High  rate  GPS  clock  estimates  applied  to  remove  SA 


7.2  Research  Contributions 


This  dissertation  provides  an  extended  evaluation  of  real-time  satellite  OD 
using  GPS  and  offers  several  suggestions  to  improve  the  state-of-the-art.  The 
affects  on  real-time  OD  of  using  precise  GPS  ephemerides  rather  than  those 
broadcast,  using  dual-frequency  ionospheric  corrections  instead  of  no  corrections, 
using  SA  free  GPS  measurements,  and  using  various  gravity  models  and  truncated 
gravity  models  are  characterized  in  terms  of  accuracy/precision  and 
computational  burden.  Also  in  this  dissertation,  GAAF  is  implemented  in 
propagation  and  OD  scenarios,  zero  bias  DRVID  is  developed  and  implemented, 
DMC  is  implemented  and  extended  to  estimate  accelerations  in  RIC  coordinates 
and  GMC  is  employed  to  reduce  the  burden  associated  with  tuning  DMC. 

7.3  Recommendations  for  Future  Work 

7.3.1  GAAF  Coefficient  Estimation 

Sophisticated  OD  filters  can  simultaneously  estimate  spherical  harmonic 
gravitational  coefficients  along  with  the  satellite  state.  It  is  recommended  that  the 
capability  of  estimating  GAAF  coefficients  in  an  OD  scheme  be  investigated. 

The  advantage  of  estimating  GAAF  coefficients  over  spherical  harmonic 
gravitational  coefficients  is  that  in  GAAF,  information  from  each  measurement 
will  only  affect  those  coefficients  whose  latitude  and  longitude  are  near  where  the 


measurement  is  taken.  In  estimating  spherical  harmonic  gravitational  coefficients 
the  processing  of  each  measurement  affects  all  coefficients  used  in  the  expansion. 

7.3.2  GMC  RIC  Extension 

GMC  has  only  been  applied  to  XYZ  DMC.  An  extension  of  GMC  to  RIC 
DMC  is  recommended  where  the  RIC  x  s  and  as  (xR,x[,Tc,aR,aI,ac)  are 

optimized.  In  this  case,  instead  of  optimizing  quadruplets,  GMC  would  optimize 
octuples.  The  Riemann  summations  would  stay  the  same  as  the  XYZ  GMC 
implementation,  only  the  x  s  and  a  s  would  be  different  for  each  of  the  three 
position  summations.  Additional  reproduction,  crossover  and  mutation 
procedures  would  also  need  to  be  added  to  take  into  account  the  extra  four 
members  included  in  each  individual. 

7.3.3  Automated  DMC  Tuning  With  the  Genetic  Algorithm 

Another  promising  application  of  a  GA  is  tuning  of  DMC  parameters 
outside  of  GEODE.  Here  RMS  and  RSS  statistics  would  be  used  to  gauge 
GEODE’s  performance.  GEODE  would  be  run  using  all  individuals  from  a  given 
generation  and  the  RMS  and  RSS  statistics  would  be  used  to  compute  the  merit  of 
each  individual.  Again,  each  individual  is  made  up  of  a  different  set  of  x  s  and 
a  s.  Then  the  GA  procedures  of  reproduction,  crossover  and  mutation  would  be 
performed  based  on  the  merit  of  each  individual  and  the  next  generation  would  be 
used  in  GEODE  again.  This  process  would  be  repeated  until  the  GA  converges 


on  an  optimal  set  of  x  s  and  a  s.  Figure  7.1  shows  a  diagram  of  the  proposed 
application  of  the  GA  to  DMC  tuning. 


Figure  7.1  -  GA  for  DMC  Tuning 


7.3.4  GPS  Point  Solution  Measurement 

There  are  several  satellite  missions  (QuikScat  and  Picasso)  where  GPS 
receivers  are  flown  but  raw  GPS  pseudorange  and  phase  measurements  are  not 
available  in  real-time  and  in  the  case  of  QuikScat  not  available  for  post¬ 
processing.  In  both  of  these  cases  more  accurate  satellite  position  knowledge  is 
required  than  can  be  determined  by  the  GPS  receiver’s  point  solution.  Therefore, 
a  modification  to  GEODE  is  suggested  where  GEODE  would  use  the  GPS 


receiver’s  point  solution  (both  position  and  velocity)  as  measurements  for  the 
filter.  The  modifications  to  GEODE  would  include  removing  the  GPS  receiver 
clock  bias  and  clock  bias  drift  rate  terms  from  the  state,  reformulation  of  the  state 
transition  matrix  and  observation  state  relationship,  and  declaration  of  additional 
variable  to  support  the  input  of  the  receiver’s  point  solution. 

7.4  Conclusion 

A  new  era  in  real-time  satellite  orbit  determination  began  when  SA  was 
turned  off  on  2  May  2000.  Meter  and  sub-meter  3D  position  accuracy  is  now 
available  in  real-time  on  board  a  satellite.  GSFC’s  GEODE  is  a  powerful  real¬ 
time  satellite  OD  software  suite  that  goes  a  long  way  to  make  meter  level  real¬ 
time  OD  possible.  Part  of  GEODE’s  power  is  realized  through  GSFC’s  provision 
of  GEODE  source  code.  GEODE’s  modular  nature  and  the  relative  ease  of 
modification  make  it  flexible  and  highly  usable  software. 

This  dissertation  offers  several  suggested  modifications  to  GEODE  to 
increase  its  accuracy/precision  and  autonomy  and  to  decrease  its  computational 
burden.  Accuracy/precision  improvements  are  demonstrated  with  the  application 
of  DRVTD,  DMC  and  GMC  but  each  increases  GEODE’s  computational  burden. 
GMC  improves  GEODE’s  autonomy  by  allowing  the  on  board  filter  to  adapt  to 
changes  in  the  dynamic  environment  but  imposes  a  significant  computational 
burden.  While  GAAF  reduces  computational  burden,  it  does  not  improve  or 
degrade  GEODE’s  accuracy/precision.  Each  suggested  improvement  provides  its 


own  advantages  and  disadvantages  and  therefore,  the  implementation  of  each 
improvement  must  be  balanced  against  mission  requirements  and  hardware 
capabilities. 
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